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ABSTRACT 


The goal of this study is to investigate the nonlinear dynamies of two ships in 
elose proximity towing. The sway and yaw dynamics of both the leading and the trailing 
ships were included in the formulation. Previous studies were restricted to a linear 
analysis, which can accurately predict the regions of stability and instability for the 
system. The mechanism of loss of stability can be assessed with a systematic nonlinear 
analysis. The analysis is based on Taylor series expansions of the equations of motion up 
to third order terms. It is shown that the primary loss of stability occurs in the form of 
Hopf bifurcations to periodic solutions. A nonlinear stability coefficient was calculated 
which allows characterization of the stability properties of the resulting limit cycles. The 
results indicate the effects of ship separation and towing tension on motion stability. 
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I. INTRODUCTION 


One of the main eoneems in both naval arehiteeture and maneuvering eontrol is 
the stability of the ship towing. In partieular, the elose proximity ship-towing problem 
has many unexplored areas of investigation ranging from building a mathematieal model 
of the two towed ships, analyzing the direetional stability of the eombined system, loeal 
stability of the maneuvering, nonlinear dynamies of the system and elastieity effeets in 
the towline. 

The close proximity ship towing is under current study by Office of Naval 
Research. Traditionally, ship towing for the military ships and the others are done with 
long towlines. The Office of Naval Research is studying the close proximity towing for 
its innovative Hull Programs, such as small waterplane area twin hull (SWATH) ships 
and their variations (such as the SLICE hull). During the SEA EANCE project, 2000, at 
the Naval Postgraduate School (Total Ship Systems Engineering Program) it was shown 
that one of the main benefits of close proximity towing in military applications is the 
ability to separate a combatant ship from a main part of its payload. 

In this thesis the focus of our efforts is on the nonlinear dynamics during the loss 
of directional stability in close proximity ship towing. Previous studies on the directional 
stability of ship towing were performed in 1964 by Abkowitz who developed the 
characteristic equation for single body towing, and by Bernitsas and others who 
developed the criteria for stability of Abkowitz 4'* order characteristic equation. The 
underlying assumption in those studies was; the leading ship is on a constant straight-line 
motion of a point mass. This is a valid approximation for long towlines while it may not 
capture the combined dynamics during towing in close proximity. Therefore, it is 
possible that the existing criteria of directional stability may be inadequate. The 
directional stability for the coupled system where both the maneuvering dynamics of 
leading and trailing ships are considered was studied by GOKCE in 2002. As a result, the 
complete system has now a dimensionality of 8 and the characteristic order is 8* order 
instead of 4*. With the inclusion of the dynamics for the leading ship, it is necessary to 
incorporate path control into the system. Eor the demonstration purposes, we include a 
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rudder control law for the leading ship only. This is done in terms of full state feedback in 
terms of the sway and yaw velocities, heading and lateral deviation of the leading ship 
from its desired track. Control law responsiveness is quantified via its time constant, 
small values of the time constant designate a responsive controller while larger values 
designate a more sluggish design. 

In this study as a leading ship we will use the SLICE vessel, which is 105 ft, and 
180 tons, and as a trailing ship the SWATH ship Kaimalino which is 89 ft, 47 ft beam, 
10-12 knots speed and 217 tons. Typical speeds for the combined vessels in the tow 
configuration are up to approximately 15 knots. The Slice is a high-speed variant of the 
Swath technology with 4 underwater hulls instead of two. Attached to each hull is a strut 
that extends up to support the main body. Figure 1. shows the profile view of the Slice.[2] 



The SSP (Semi Submersible Stable Platform) Kaimalino was the world’s first 
high performance open ocean Swath ship. It consists of two parallel torpedo-like hulls. 
Attached to the hulls are two streamlined struts. The struts extend above the water surface 
and support the main body. The Kaimalino also has stabilizing fins attached near the aft 
end of each hull. Figure 2 shows a profile of the SSP Kaimalino. [2] 
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Figure 2. The Profile of the SSP KAIMALINO [From; [2]] 

The general approach adopted in this thesis is as follows: 

1. First we will use the maneuvering equations of motion in the horizontal plane for 
both leading and trailing ship. 

2. Coupling between the two ships will be provided through the towline. 
Hydrodynamic coupling arising from radiation and diffraction effects will be 
neglected. This can be easily incorporated into the analysis once the effect of such 
hydrodynamic coupling on the hydrodynamic coefficients is established. 

3. The system will be linearized in order to study its stability properties. At the end 
of this part of study identification of regions of stability is expected. 

4. A nonlinear analysis will be performed for identifying mechanism of loss of 
stability. The response of the system will be analyzed with the help of both the 
linear and nonlinear studies. 

The hydrodynamic coefficients of both ships used in this study are provided from 
reference [2]. 
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II. PROBLEM FORMULATION AND STABILITY ANALYSIS 


A. EQUATIONS OF MOTION 

In previous studies, several different values for tension and for length of the 
connection between the two ships were tried for the stability analysis of the ship towing. 
[Gokce, March 2002] For every instant, a zero eigenvalue was found. This means the 
system cannot reach to stable condition. This zero eigenvalue was attributed to the lack of 
directional stabilizing effects on the leading ship. When the signs of the towline restoring 
forces and the moments are investigated, it can be seen that towing force has a 
destabilizing effect on the towing ship. On the other hand it has a stabilizing effect on the 
towed ship. The ships cannot show directional stability in the horizontal plane, therefore 
a rudder control law must be applied to the equations of motion. 

The mathematical model of the surge, sway and yaw maneuvering equations of 
motion in dimensionless form is shown below. The subscript 1 refers to towing (leading) 
ship and the subscript 2 refers to towed (trailing) ship. In these equations u is the surge 
velocity, v is the sway velocity, R is the resistance of the vessel moving through body of 
water, and T is the tension in the connection (rope, cable, or some other mechanism) 
between the two towing ships. The rudder angle is denoted by 6, the heading angle is \|/, 
and Y is the geometric angle depicted in Figure 3.The equations of motion are: 


(^2 - r.2)V2 - “ T sin(^^2 + 7) (1) 

{Iz2 ~ ~ ^v2^2 ~ ^v2^2 ^r2^2 ~ p2 Sill((^2 7) (2) 

(«! -4)Vj - T./j = -m,r, + T^jVj + + T sin(^^j + /) + YgS (3) 

(4i - )h - NJi = KiY + - TXpi sin(^^i + /) + NgS (4) 

where 

^ - 4ti (5) 
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To study stability of this system, the equations must be linearized. We are dealing 
with small heading angles, so the small angle assumption ean be made. 

sin 1^1,2 =1^1,2 (6) 

eos^//j2=l (^) 

Another assumption we ean make is for the veloeities of the ships: 

Mj = = 1 (8) 

Therefore, the linearized set of equations of motion 

(^2 - Y,2 >2 - Y.^r^ = -^2^2 + Y,2V2 + Y,2^2 “^(1^2+ j) (9) 

{ 1^2 r 2)^2 ~^ v 2^2 ~^ v 2^2 ^ r 2^2 ~ Yx ^ 2^2 ( 10 ) 

(mi -4)vi -T,iri =-miri +T,iVi +T,iri +r(^^i +y) + T^^ (11) 

(4i =iV.,Vi + A^iri -7x^l(^^l +y) +(12) 

The kinematies and the geometry relations of the system are deseribed by 
equations (13) and (14). 

jj = Ml sin^/Zj + Vj oos(^i (13) 

j>2 = Mj sin^/Zj+^2 cos(^2 (14) 

If we define 

T = T2-Ti (15) 

then 

T = T2-Ti (16) 
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These geometrie relations are explained in Figure 3. 



When the equation (13) and equation (14) are substituted into equation (16), we 


get 


y = U 2 sin y/^ + V 2 eos y /2 - sin y/^ - eos y/^ 
From the Figure 3 it ean be seen that; 

T2 + ^p2 sin y /2 - (ti - Xpi sin Wy ) 


(17) 


sin 7 : 


/ 


(18) 


whieh ean be rewritten as 


sin7 = ^ + y(x^2sin^^2+x^isin0) (19) 

Using the small angle assumption onee more we get 

Ti=i^i+U (20) 
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y^ = y/^+v^ 


( 21 ) 


y 2 +x ^y/^-(y,-x ,y/,) 

sin y =- - - - - 

/ 

The next step is to eonvert these equations into a matrix form. 


'^1 ~ ^2w^2 ^Ivr^l ^2v iW 2 T) (23) 

^2 “ ^2rv^2 ^2rr^2 ^2r ^2 T) (24) 

n = + 7 )+ C,^S (25) 

= An,Vl+Arrri+BiXWl+7)+Ci,S (26) 


The eoeffieients in equations (24), (25), (25) and (26) are found after the 
following mathematieal steps: 


4„ = [(l.2-NyY„+(N„Yy]^[(m,-Y„ll,,-Ny-N„Yj (27) 

4.. =[fc -A',2)-4.2r2] (28) 

4.=[-(/,2 -A',2)?’-(r2?’42)J-[("<2 -A',2)-y.2r2l ( 29 ) 

A„=[{Y.2N„)yN,M 2-Y,2M(m,_-Y,Jl^,-N,.2)-N„Y„] ( 30 ) 

4.. =[(r2 -”<2)4.2 +42('«2 -r.,2)l-[(/«2 -r2X/-2 -42)-42r2] (3 1) 

4. = l-r42 -r2)J-[(”<2 -r2)(2,2 -42)-42r2] (32) 

4» =[(32, -4,Ir, +(4,r,)]-[(<”■ -4,X4, -4,)-4,r,] (33) 

4.. . = [(4 -»<,)(3„ -4, )+(4,r, )]-[(”<■ -n,X4, -4,)-4,i'„l (34) 

4. = 1(32, - 4, )?■ - (r,n .-,2 )f [(”<■ - Y., ti,, - 4,) - 4,r ■ 1 ( 35 ) 
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=[fc -»>,)»., +iV„K -nJl-K"-, -n,X4 (37) 

B„ =[tN„+Tx,M, -i’,)l-[(»>, -i’lX/., (38) 

C„ =[(/., -A'„fc + A',3’„]-[(m, (39) 

C„ = fciV„ +A',(m, -F„)]-[(m, -F„)(/„ -Af„)-iV„F„] (40) 


The hydrodynamic coefficients are taken from reference [2] with the exception of 
Yg and Ng. There is no data available for the rudder design of the SLICE, so we 
assumed a ship turning radius of three ship lengths under fifteen degrees of rudder, and 
calculated the necessary values of the rudder hydrodynamic coefficients to achieve that 
turning radius. 

The control gains k^,k^, k^and k^ are calculated by standard pole placement 
techniques by utilizing the Matlab program shown in the appendices. 

The system matrix of 8x8 for the stability analysis is as follows: 


1^2 


1^2 

^2 







= M 


T2 


T2 

Ti 


Ti 

Wi 


Wi 

_¥i\ 




where 


(41) 
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M= 


0 0 


A. 



X 

Av + 


l 

i 


/ 

1 





^2r^p\ 

1 

/ 


A.+ ^ 

1 

A. 



^iv^p2 

A.x„ 

/ 

l 


1 

l 



r k 

B\rX p2 


i 

l 


1 

/ 

0 

0 


1 

0 

0 

0 


0 

1 

0 

0 


0 

0 

0 

0 


0 

0 






(42) 


Once the matrices are formed, we ean predict the stability of the system by an 
eigenvalue analysis. This analysis is outlined in the next seetion. 


B, DEGREE OF STABILITY 

The eritical eigenvalue of the system is the one that has largest real part. This 
largest real part is the degree of the stability of the system. If we have a negative real part 
it means we have stable system. On the other hand the eigenvalue with a positive real part 
indicates that the system is unstable. 

C. REGIONS OF STABILITY 

The results ean be summarized in a single graph designating regions of stability 
and instability. This graph is shown in Figure 4, whieh was produced by utilizing the 
Matlab eode shown in Appendix A. The graph shows the eritieal value of the eontrol time 
eonstant for stability versus the towline length, and is parameterized by the tension in the 
towline. Combinations of (Tc,L,T) below the eorresponding curve will yield stable 
response, while combinations that are located above the curve will result in system 
instability. 
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Figure 4. The Stability Region 
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III. NONLINEAR STABILITY ANALYSIS 


A, HOPF BIFURCATION 

In all cases of the loss of stability that we defined in the previous chapter, a pair of 
complex conjugate eigenvalues has zero real parts. This is called Hopf Bifurcation and is 
accompanied by periodic solutions (limited cycles). The Hopf Bifurcation occurs when a 
pair of complex conjugate eigenvalues crosses transversally the imaginary axis. When 
this happens, the system will deviate from a steady-state solution in an oscillatory 
manner. This deviation can be either supercritical or subcritical. For the supercritical 
situation, stable limit cycles form after straight-line stability is lost. Periodic solutions 
appear after the initial loss of stability and the resulting periodic solutions are stable. For 
the subcritical cases, the periodic solutions appear before the initial loss of stability and 
the resulting periodic solutions are unstable. In such a case, the domain of attraction of 
the equilibrium point becomes progressively smaller as the critical point is reached, and 


c 



Figure 5. Supercritical and Subcritical Poincare-Andronov-Hopf 
Bifurcations 
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in fact it shrinks to zero at the eritieal point. Therefore, the equilibrium point is unstable 
under finite disturbanees and in praetiee it may lose its stability even though the linear 
analysis is predieting it to be stable. Classifieation of the nature of these bifureations is 
the subjeet of this ehapter. In order to evaluate the stability of the limit eyeles of Hopf 
Bifureations, a nonlinear analysis is needed to be performed. 


B. FORMULATION 

In the previous ehapter, in order to get the degrees of stability and the regions of 
stability, the equations of motions were linearized. In this way, the linear stability of the 
ship towing was studied. In this ehapter the nonlinear dynamies of the system will be 


studied. 

The nonlinear forms of the equations were deseribed before: 

(«2 - ^ V 2)^2 - ^, 2^2 = -^ 2^2 + Y ^ 2^2 + ^ r 2 fi “ ^ Sin (^^2 + 7 ) (1) 

(■^Z2 ~ ^f2^^2 ~ ^v2^2 ~ ^v2^2 ^r2^2 ~ ^p2 Sin(^2 ^ (2) 

(mi - T-i )vi - T-iq = -m,r^ + + T sin(^^i + y) + Yg5 (3) 

(4i - -TXp, sin(5^i + /) + NgS (4) 

When we write these equations in the matrix form; they are as follows: 

1^2 = ^2vv^2 + + ^2v sin (l^2 + Y) ( 43 ) 

h = Arv^2 + + K sin {¥2 + 7) ( 44 ) 

+ A.A + K sin ( + 7) + (45) 

h = ArvY + ArA + Ar sin {Vi + Y) + C,,d ( 46 ) 

y = U2 sin ¥2 + Y ¥2 ~ sin ¥i ~ ^1 cos ¥\ ( 47 ) 
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The first step for the nonlinear dynamics analysis is to expand the y term in 


Taylor series up to the third order. We have: 

y^+x^smy/^-[y,-x,smy/,) 

smy = - - -- - - (18) 

by utilizing the sin\|; taylor expansion: 

1 3 

sm^ = ^ — y/ (48) 

6 


y=\ll 


yi+Xpi ¥iy¥i -Ji+Xpi ¥x-y¥\ 


(49) 


This alpha value is applied to the equations (43), (44), (45), (46). In order to 
calculate the nonlinear dynamics coefficient, we have to form a matrix of g values. By 
these we refer to the nonlinear terms in the equations of motion up to third order. These 
g (x) values are obtained from the equations (50), (51), (52), (53), (54). 


^2 = 4w^2 + Av/i + K sin(^, + ;y) (50) 

giW 

4 = Arv^l + Arr^l + K + (5 1 ) 

g2(x) 

+ AvA + Av sm(l^i + r) + Ci,(^ (52) 

gsix) 

h = + ArA + Bi,sm{iy^ + y) + C^^S (53) 
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( 54 ) 


y = sin y/^ + Vj cos y/^ - sin y/^ - Vj cos y/^ 

“V ■ ' ' ' ' 1 — 

With these g (x) values a matrix of 1x7 is formed: 


g(x) = 


gi(x) 

gii^) 

g^M 

g4(x) 

gsix) 

g6(x) 

giix) 


(55) 


Next step is to apply the Taylor series expansion for the sin(\|/ +y) part of 

equations (50), (51), (52), (53) and keep the nonlinear terms up to third order. The terms 
gg(x) and g 7 (x) are equal to zero. In order to aecount for rudder angle saturation we 
substituted the discontinuous saturation function with a smooth function based on the 
hyperbolie tangent of the rudder angle. Sinee we are dealing with nonlinear analysis, the 
nonlinear terms up to third order are kept and the following coefficients are obtained. 


1 3 .,,3 1 2.,, , ^ 2 


2L 


1 2 1 2 1 1 2 

2 y 2^^ 1 y iV 2 n, t3 kl T 2 ~^p2y^2 ^ ^2 Tl 1^2 


2U 


2U 


L 


2V 


-^Kl^p2¥l¥2 -^^pxK2¥l¥l+:^xl,Xp2¥l¥l +^TlT2^p2?^2 -^TiT2^P2?^ 

1 1 3 1 1 3 

+ 7rTiT2^pi?/^i - 77rTiT2^p2?/^2 + 7rTi^pi^p2?/^i?^2 - 77rTi^pi^;,2?^i ¥2 

L oL L oL 

~ y\^p\^p2¥\¥2 2,^1] yy^y2^pV^p2¥\ ¥2 y2^pV^p2}¥W2 


y^y2^px^p2¥l¥l-^y2^px^p2¥,¥2-^^pi¥,¥l+:^^px¥l¥l-^y2^p2¥l 
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+ ^yx^p2¥l + jTyiy2¥2-i^^lx¥i¥2-jjy2¥l + jjy,¥l-i^xl^¥l 

1 3 2 1 2 1 1 3 1 

^-^’^pl^p2¥l¥2 —^XpxXp2¥l¥2 —^y2XpxXp2¥l¥2+-^y2Xpx¥l¥2+-^yiXpl¥,¥2 

OL L L oL L 


2 .,,2 1 2 .,,2 1 2 .,,3 1 2 


g£2 ¥2 2jJ> y2^p2¥2 + 2^3 y2^p2¥2 y2^p2¥2 


-^y\^px¥,^^yl^p,¥l^^yA2¥l-^yl^p2¥2^^yUp2¥\ 
-^fi^pi¥i +:^ylxpi¥l-^y2Xp,¥f -^v4-^Xp2¥l +\xpi¥i-^^ p^¥\ 

^ 3 „ J , 1 3 , 1 2 .,, 3„,2 1 33 


^2 = ^2. {-:^yi^pv^p2¥\¥l ^^^\v^p2¥l¥\ +^y' -^y2¥2 -i^yl¥2 

JOL VIL oL IL IL 


1 .,,2 1 2 .,,2 


+—yx¥2 - i^^px¥, ¥2 + -T7Ty2^pi¥, ¥2 - -vTy2Xpi¥i¥2 + ^yiy2¥ 


+ -^^pi^p2¥l¥l - ^XpiXp,Y,¥l - ^¥l ~^y"v^2 + ^yi^px^p2¥i¥2 
--^yiXp,Xp,wf¥2-^xl,Xp,w^i/^2-^Xpyp2¥i¥2+Y^Xpixl2¥i¥2 


+ :^yyp2¥2 + j2yi^p2¥l -i^yl^pi¥i + i^ylx^,Yl + ■^ylx,^¥l 


-^^p,¥i¥l -I^y2^l2¥l - -^yiXp,¥l¥2 + ^yi^pi¥,¥2 - -^^pi¥l 

IL IL oL L oL 
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1 1 1 3 1 1 3 

~r^\~r^piV^2 y2^piV^\ ~TTy2^p\^p2y^2~7T2 y2^p\^p2y^iV^2 

L L oL L oL 


+ 


1 


y 2 ^ pi^ p 2 y^iy^ 2 


1 

3^ 


y 2 ^ p\^ p 2 y^ \ 2 


1 


y'y2 + 


1 


J 1 J 2 - 


1 



3 

2 


6^3 -^2 2£2 ^ p2W 2 ■*■ jJ, y\y2^p2W2 


1 3 1 

-^yxy2^p2¥2 + -^yxy2^px¥i 


1 


^\x¥ 


3 

1 


1 


21 


^px¥i + 


121' 


-yl^px¥l + 


IL 


3 


31 




3 

2 


+ -^^pi¥\¥l -^y2^p2¥l 


Z2 - -^72V: 


1 .2 1 


1 


+ — J2?/^r + -^y2^px^p2¥x¥\ - -J2yi^pl¥l - ^^p2¥l¥2 + p2¥l¥l 


^^yx^px¥l - -^^\2¥x¥l + ^ Jl J2?^l + J^^pl^p2?/^lV2 - -^yl^px¥x 


+ - i^y - -ATyi^pi^p2¥i¥2 + 7^^li^p2¥i¥2 


2L 


2U 6L 


\2L 


^^r^yi^pi^p2¥l¥l + ^Xp,Xp2¥l¥l-^Xp,Xp2¥W2-^yxy2^p2¥l 

5 oL oL L oL 

"*" ~f3 y \y 2^ p2¥ 2 "*" r t 2 y 2^ p2¥ \¥ 2 "*" ~ 7 T¥ y 2^ p\^ p2¥ \¥ 2\^ p2¥ 2 "*" T^TTF 3 " 1 ^ p2¥ 2 
L oL oL iZL 


1 


jjYyf^px¥i + Yh^y'^P^v^' + :^yxKx¥x 
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1111 1 3 

-^yiXpiWiWi+^yxyi^piWx +^yxXp2¥x¥2 -^y2XpiXp2WxW2 --^yxy2^px¥x 


+■ 


Is 


-yx^pX^p2¥xV2 


361 


3 y2^pX^p:^X ¥2 


2L 


2 3 

-Xpx¥x 


3L 


Xpx¥x + 


2U 


■Ji72 


^^lx¥l -^^P2¥l-^yly2-\¥l-^yx^p2¥x¥l -^^\x^p2¥l¥2 

OL OL IL O OL IL 


1 1C 

j^yxXpxXp2¥l¥2 -s^x^,xl,y/,Yl}--^{klKwl + klk^y,Yx +K^yx¥x 


+k k^ry + k k\y + 2k^k^k v^r,y^ + kyy + 2k k^k^v.r^ + 2k k^k v^y.y/, 


+'^k k^k ry,Y, + klk^r^ + k^^k vfy^ + kXr^ + kXvX + k^kXy^ +^kXi 


2 „ 2 , 


2 ,. ,.2 , i 2 i „ 2 , 


1 


- 3.„3 


+\ks+\ ey+\kiy]) 


&= «S-XXyl~yy, v.V ,-XwiyXwl 

+ ^y2^px^p2¥x¥l-^y2^px¥l -^^p2¥X2+:^^p2¥XX^yx^px¥x 


1 2 2 1 


2L 


2 X¥x¥;+j2yxy2¥x + kwvy^pxyy¥2 


1 

121" 


2V 


-y2^px¥x 


+ ■ 


121 - 


■y2Xpx¥x 


' y2^p2¥2 22^3 y2^p2¥2 ^ j2 y2^px¥X ^ ^3 y2^pl}!^2 ^ 7^3 yX^p2¥2 


\2U 


2U 


2U 


2U 


1 


2L 


? 1 7 ^ 1 9 

TyxXp2¥2+ . 


^2^3 yx^p2¥2-i^y:^px¥x + i^yxXpx¥x+^yxX¥x 


^yx^px^p2¥x¥Xi^Xxp2¥l¥X^yx^px^p2¥l¥X^Xpx^p2¥Xl 
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1 2 1 3 1 1 3 

j^2 ^p\^pi^ 1 2 y\y2^piW 1 ■*■ ~jy y\y 2 ^p2y ^2 ^^2 y 2 ^p2V\V 2 

1 j 1 1 1 1 1 j 

-^yi^piWWi ^^yxyi^pxWx ^^yx^^iwy^i-^yy^px^piWxWi —^yiyy^pm 
^^yi^px^p2Vx¥2-^^y2^px^p2v\v\-^^\xV\-^^px¥l 


-■^^;, 2?^2 -^ 71^2 ^ ^ Ix ^ p2 ¥^ x¥2 


1 1C 

-J^yx^px^p2¥l¥2 -l^^pX^l2¥x¥l}--^{KKrx¥x+K^yyx¥x+K^yx¥x 


+k kW¥x + KKyy + ykrkyWxyx + kyy + Ik k^k^v.r^ + Ik k^k v^y.y/, 


+^K^rkyrxyx¥x + Kkry + + KK^X + +^K¥x 


2 ,. ,.2 , i 2 i „ 2 , 


3 V 1 3^1 ^ yJ-Xi 


In order to get g^, a Taylor series expansion of sine and eosine is used. 

1 3 

smY = Y~~¥ 

1 1 2 
cosy-^~ — ¥ 

g5=¥2-¥x-\(¥l-¥l ) + ^2 -^1 - \y¥l - Wl ) 

The next step is to transform eoordinates to normal form. This is obtained by the 
using matrix of eigenveetors [m] at the eritieal point as the eoordinate transformation 
matrix. The new set of coordinates is called [z], so in terms of the original set we have 
[x]=[m][z]. The only variables of [z] that need to be retained in the analysis of this 
problem are first two, the remaining variables are expressed in terms of first two 


20 



variables through higher order relationships and do not affeet the third order expansions 
that we utilized. 




Vj =m2l2l+«22^2 


h="^3A+«32^2 


=m4iZj+m42Z2 


¥2 ^51^1 ^52^2 


V2="^61^1+"^62^2 


Tj =m2iZj+m72Z2 


y2=Vx+V2 


When the above terms are substituted into the equation of motion, we get; 


1^3 3 33 33 ', 1 / 32,2 


^2v{^j2^^U ^81 ^5\^p2 ^11-^pl) 22^2 l-^51-^p2 ^41^51 ^11^51-^pl) 


6U 


12 1 1 




-mi>5iX^lX^2) - JT(f^5lf^&l^p2 - «^41«^51^p2 “ 11 5 1 ^ 8 1 

-'«ll«4l"^51^pl +"^ll"^5Vpl^p2)-^K"l"^81 -"^4X1 +'”lXl^pl)+^Kl"^5l"^81^;,2 
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^\2 ^2v{g^3 (^42 ^82 ^52^p2 ^12^pl) 22^2 (^52^p2 ^12^52^pl) g ^52 ^j^^52^p2 


^12-^pl '*’ 22^3 ^^42^82 ^42^82 ^52^&2^p2 '*’ ^42^52^p2 ^^12^82-^pl ^l2^S2^pl 

-mi2«5V;,l^j2 -"^f2"^52^Jl^p2)-^("^82"^52^p2 
+ 2m42"^52 “ "^42"^52"^82 + "^12"^52"^82^;,1 “ "^12"^52"^42^;,1 + "^12"^5\^pl^;,2) 
-^(^52^82 - 'W42»252 +»2l2»*52^pl) + ^(^42»«52^82^p2 + »2l2»242^82^pl “»«12»«52»*82^;,l^p2 


+m|2«52«42^pl^p2) 


Aa -^ 2 v { 22^3 ^^41^42 ^81^82 ^5\^52^p2 ^ 11 ^ 12 -^pl '*’^42^81 ^41^82 ^5\^S2^p2 

+m^^ml,x]^ - 2mi>82xji - + ml.m^^x], - m.^ml^x^.x]^ 


mi>52x" X„ 2 ) - ^(^«^51«^52^A - »241«^81^82 + m + ttl ,,m ^,X 




U 2 


5r'^52*^p2 ''Mr'^8r'^82 ' ''Mr'M2''^81 ' ''^5r'^52”^81*^p2 


-«4l"^51 V^2 +Mlf^2f^Al +«ll"%l Vpl +«4l"^42"^51^p2 +"^l"^41 V;,l -^1^2^]^ 
-"241^51 Vp2 -»242»251^1^p2 “»2l2'”41»%l^pl "l»241»%2^pl -»2ll»242»*81 +»2ll»252»281^pl^p2 
+ ^12^51^81X^1X^2 + mu^ 5 l^i 2 XplXp^ - ^11^41^52X^1X^2 - mX 2 ^ 4 im 5 lXplXp 2 

-»2ll«^42^5lXplX^2 +«^11«^51«^52 Vp 2 +«^11^12^51^pl^p2)-^^5>52 -^K>12^;,1 


+ml,m^^-m,^ml,+m,^ml,Xp,)-jj{2m,,m,^m^,Xj,^+ml,m^^Xp^-2m,,m,,m,^Xp^ 


-"^42"^51^p2 + «41«42"^51 “ «4l"^52"^81^p2 “ «4l"^52"^81 “ "^41«5l"^82 “ «42"^51«81 
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+^11^52^81X^1 + ^12^51^81X^1 + ^11^51^82X^1 - ^11^41^52X^1 - 


Z 2 1 

-^1 im42/«5lX^l + 2^1 im5im52X^lX^2 + «^12^51^pl^p2 +l«^12^51^pl ) - - Ki^52^81 


-m4im5im52 +miim5im52X^i)-^(m4>52 + m4>52X^2 + <i"^52^pi) 


^14 ^^'’^'2V’ ^81^82 ^5\^52^p2 ^\2^ p\ ^ 41^&2 ^42^81 


-«52"^81^?2 + «41«5V?2 - "^42"^51^;,2 - "^n"^4Vpl + «f2"^41^?i - "^u"^5Vpl^p2 


1 3 

~^12^5lX„iX 2 ) 7T^~Z^51^52^p2 '*’ ^52^81-^n2 '*’ ^^51^52^82-^n2 ~ 


U 2 


-2^42^51^52X^2 - 

+ ^ 11 ^ 52 ^ 82 X ^1 + ^ 12 ^ 51 ^ 82 X ^1 - tU .^tU ,^X - ^ 11 ^ 42 ^ 52 X ^1 - ^ 12 ^ 42 ^ 51 X ^1 


+ miim52X^lX^2 - 2»2i2/«51»252^pl^p2 + 112 52 ^pl + H 5 ) “ “ »« 5 152 


- ^(2m5im52^p2 + «^51«^52»282 - »«42»251«^52 + «^12«^51»«52^pl) - 

+ »«52^81 -^41«^52 +«^11^52^pl) + ^(«^42«^81«^82 -^41^42^82 -«^51^52«^82^p2 

+ "^42«51"^52XJ2 - "^4l"^42"^52^p2 “ "^12"^4l"^42^;,l + "^1 l"^12"^42^pl + "^4l"^52"^82^p2 
+ ^42^52^81X^2 + »242^51»282^;,2 + »2l2^41^82^pl + »*12»242^81^pl + l»242^82^pl 

-^ 12 ^ 52 ^ 82 X^ 1 X ^2 - »2ll»252»282^pl^p2 - »2l2»251»282^pl^p2 + 12 41 52 ^pl^p 2 


+ miim52m42X^lX^2 + ^n^AZ^iX^pl^p2 “ ^n^Sl^SZ^pX^lz “ "^ll"^12"^52^Jl^p2) 
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Following the above background calculations we are now ready to summarize the 
nonlinear analysis for our system. 


C. THE CONTROL LAW 

A typical full state feedback control law will be utilized for the nonlinear analysis 
again. In the below equation So is the feedback rudder angle for small deviations. 

= k.W,^k.v.^kr+k.y, 

In order to capture the effect of rudder saturation, the control law is modeled as; 




C \ 

do 


\0 sat J 


(57) 


The hyperbolic tangent function is used instead of hard saturation function 
because of its analyticity properties. 

Utilizing the taylor expansion again; 


1 

smy/ = y/- — \f/ 


3 


S - do- 




Si 


(58) 


In the Matlab (Appendix A and B) code for the nonlinear analysis saturation limit 
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on the 5 is set to 0.4. 



D, INTEGRAL AVERAGING 

When the equations of motion are written in the compaet form, 

X = f{x), X = (59) 

The system can be written in the form 

X = AX + g{x) (60) 


A is the Jacobian matrix of f (x), and g (x) contains all nonlinear term of 
equations (50), (51), (52), (53), (54). 

If we introduce the transformation matrix (T) of eigenvectors of A 


T-' = K] 

i,j = 1,2,3,4,5,6,7,8 

(61) 

T" = [n,] 

i,j = 1,2,3,4,5,6,7,8 

(62) 

linear change of coordinates 



x = Tz 

z = T~^x 

(63) 


will transform system (X) into its normal form 

z = T~'ATz + T-^g(Tz) (64) 

with this transformation we get 

O-wO 0 0 000 

-wOO 0 0 000 

0 OCjOOOOO 
0 00 C, 0000 

(65) 

0000 C3 000 

0 000 OC 4 OO 

0 000 OOC 5 O 

0 000 000 Cg 
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The expressions for g are given by: 


/jlZj + li2^2 ^2 hA^\^2 

^21^1 ^22^2 ^23^1 ^2 ^24^1^2 

g — ^3 2-^2 ^3 3-^1 ^2 ^4^1^ 2 

^41-^1 + ^42^2 ^43-^1 ^2 ^44^1^2 

h\^l “*“ ^52^2 “*“ ^2 “*“ ^54^1^2 J 

(66) 

If we substitute equation (63) into (60), we get 


z, = -wzj + Tiiz' + Tijzfzj + qjZjZj + q^Zj 

(67) 

Z2 = —WZj + TjjZj + TjjZj Zj + r23ZjZ2 + ^24^2 

(68) 

These equations ean be written as, 


Zi — —U’Z2 +(Zj, Z 2 ), 

(69) 

Z 2 — “WZj + F 2 (^1 ? ^ 2 )’ 

(70) 

Introdueing the polar eoordinates for the equations (69) and (70), 


Zj = 7? eos 0, ^ 2 - ^ ^ 

we get 

(71) 

R = F^{R,0) eos G + F2{R,0) sin 0 

(72) 

RO = wR + F2{R,0 )co?,G - F^{R,0)smO 

(73) 

Equation (73) gives 


R = P(0)R^ 

(74) 
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The function P{d) is a 2;r periodic function in the angular coordinate 0. In 
order to get the equation with constant coefficients, equation (74) is averaged over one 
cycle in 6. This way we get 


R = KR^ 


and 


(75) 


K = 


— Vp{e).d0 
2^ Jo 


(76) 


After making the integral calculations we get the simplified form of the equation (76) 


^ = -[ 3 fi+f 3 +''22+3^24] 


(77) 


The indicated coefficient are given by, 


hi ~ ^11^11 ^12^21 ^ 134 i ^ 14^41 ^IS^l 

63 ~ ^1 i 44 ^ 12^24 ^13^4 ^14^4 ^1544 

p2 ~ ^ 2 i 43 ^22^23 ^23^3 ^2443 ^25^3 

p 4 ~ ^ 2^12 ^ 22^22 ^ 23^2 ^2442 ^ 25^2 


(78) 


The coefficients are the elements of the inverse of transformation matrix T, 

which was defined as the matrix of eigenvectors of the linearized matrix A. 

Results in terms of the nonlinear stability coefficient K are shown in Figure 6. It 
can be seen that the short towline lengths result in supercritical bifurcations, while larger 
lengths result in subcritical bifurcations. The implications of this on system stability 
under finite disturbances were discussed previously. It should be pointed out that the 
transition from supercritical to subcritical bifurcations as towline length is increased is 
dependent on the hydrodynamic characteristics of the two ships and is not necessarily 
true for all ships. The procedure that we have established in this work should be followed 
in order to assess the type of corresponding bifurcation. 
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Kcoef 



L 


Figure 6. Nonlinear Dynamies Coeffieient vs. Towline Length 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

This thesis presented a eomprehensive study of the stability region properties and 
the nonlinear dynamies response of the two-body ship-towing problem. The results of 
this study ean help in establishing rational guidelines, whieh should be followed in order 
to ensure stability and safety of elose proximity ship towing operations. The main results 
from this study ean be summarized as follows: 

1. The traditional towing stability model eannot prediet the dynamies of towed ships 
in elose proximity. One of the reasons of this result is exeluding the dynamies of 
the leading ship in the overall analysis. 

2. Path stability depends on elassieal parameters sueh as tension and towline length, 
but also on additional terms sueh as the leading ship’s eontrol time eonstant. 

3. For very large values of the towline length, the effeet of the leading ship eontrol 
beeomes less pronouneed. Therefore, previously reported results in the literature 
ean be eonsidered as a large towline length of the two-body model developed in 
this work. 

4. The primary meehanism of loss of stability is in terms of Hopf Bifureations to 
periodie solutions. The bifureation ean be elassified as either supereritieal or 
suberitieal in terms of a eompatible nonlinear stability eoeffieient. 

5. Hopf Bifureations tend to be supereritieal for small values of towline length, 
while they tend to turn into suberitieal as the towline length is inereased. 

B, RECOMMENDATIONS 

Reeommendations for eontinuing studies are the following: 

1. Analyze the properties of stability and bifureations ineluding the towline elastieity 
into the model. 
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2. Analyze the stability properties of different vessels and investigate the oeeurrenee 
of statie bifureations. 
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APPENDIX A. MATLAB CODE TO FIND THE REGION OF STABILITY 


% Stability graph - TC vs. L - constant T 
"6 

% Parameters: 

% T = Nondimensional tension 

% L = Towline length 

% TC = Control law time constant 

% xpl = Attachment point (leading ship - positive forward 

of amidships) 

% xp2 = Attachment point (trailing ship - positive aft of 

amidships) 

"6 

% Constants 


ul 

= 

1; 

u2 

= 

1; 

ml 

= 

0.018078; 

m2 

= 

0.018; 

Izl 

= 

0.0007; 

Iz2 

= 

0.00069412; 

Yvl 

= 

-0.07893; 

Yv2 

= 

-0.1183; 

Yrl 

= 

-0.004044; 

Yr2 

= 

-0.0042; 

Nvl 

= 

-0.016428; 

Nv2 

= 

-0.0187; 

Nrl 

= 

-0.010332; 

Nr2 

= 

-0.0176; 

Yvdotl 

= 

-0.051328; 

Yvdot2 

= 

-0.0184; 

Yrdotl 

= 

0.005617; 

Yrdot2 

= 

-0.0011; 

Nvdotl 

= 

-0.001945; 

Nvdot2 

= 

-0.0008489; 

Nrdotl 

= 

-0.00564; 

Nrdot2 

= 

-0.0090; 

Ydel 

= 

0.0103; 

Ndel 

= 

-0.0051; 

indexL 

= 

0; 

index 

= 

0; 

xpl 

= 

0.5; 

xp2 

= 

0.5; 


"6 

% Enter constant T 


T = input('T = ' ) ; 

"6 

% start loop on length 
"0 

for iL =0.1:0.01:2.0; 

indexL = indexL+1; 

indexTC = 0; 

L = iL 

L_v(indexL)= L; 

A3 = 1/L; 
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B3 = xpl/L; 

C3 = xp2/L; 

D3 = -1/L; 

g, 

0 

% Loop on TC 

g, 

0 

for iTC =0.1:0.01:2.0; 
indexTC=indexTC+l; 

TC =iTC; 

bpole = -1/TC; 

TC_v(indexTC)=TC; 

g, 

0 


% Setup the matrix coefficients 

g, 

0 



A2vv 

= [((Iz2-Nrdot2)*Yv2)+(Nv2*Yrdot2)^ 

]/ 

[ ( (m2- 

Yvdot2) 

*(Iz2-Nrdot2)) 

-(Nvdot2*Yrdot2)]; 




A2vr = 

[((Yr2-(m2*u2))*(Iz2-Nrdot2))+(Nr2*Yrdot2)^ 

]/ 

[{{m2- 

Yvdot2) 

*(Iz2-Nrdot2)) 

-(Nvdot2*Yrdot2)]; 




B2v 

= [-((Iz2-Nrdot2)*T)-(Yrdot2*T*xp2)^ 

]/ 

[((m2- 

Yvdot2) 

*(Iz2-Nrdot2)) 

-(Nvdot2*Yrdot2)] ; 




A2rv 

= [(Yv2*Nvdot2)+(Nv2*(m2-Yvdot2))^ 

]/ 

[((m2- 

Yvdot2) 

*(Iz2-Nrdot2)) 

-(Nvdot2*Yrdot2)]; 




A2rr 

^ [((Yr2-(m2*u2))*Nvdot2)+(Nr2*(m2-Yvdot2))^ 

]/ 

[{{m2- 

Yvdot2) 

*(Iz2-Nrdot2)) 

-(Nvdot2*Yrdot2)]; 




B2r 

= [-(T*Nvdot2)-(T*xp2*(m2-Yvdot2))^ 

]/ 

[{{m2- 

Yvdot2) 

*(Iz2-Nrdot2)) 

-(Nvdot2*Yrdot2)]; 




Alvv 

= [((Izl-Nrdotl)*Yvl)+(Nvl*Yrdotl)’ 

]/ 

[((ml- 

Yvdot1) 

*(Izl-Nrdotl)) 

-(Nvdotl*Yrdotl)]; 




Alvr 

= [(Yrl-ml*ul)*(Izl-Nrdotl)+(Nrl*Yrdotl)’ 

]/ 

[((ml- 

Yvdotl) 

*(Izl-Nrdotl)) 

-(Nvdotl*Yrdotl)]; 




Blv 

[((Izl-Nrdotl)*T)-(Yrdotl*T*xpl); 

]/ 

[{(ml- 

Yvdot1) 

*(Izl-Nrdotl)) 

-(Nvdotl*Yrdotl)]; 




Civ 

= [((Izl-Nrdotl)*Ydel)+(Ndel*Yrdotl)’ 

]/ 

[((ml- 

Yvdot1) 

*(Izl-Nrdotl)) 

-(Nvdotl*Yrdotl)]; 




Alrv 

= [(Yvl*Nvdotl)+(Nvl*(ml-Yvdotl))’ 

]/ 

[((ml- 

Yvdot1) 

*(Izl-Nrdotl)) 

-(Nvdotl*Yrdotl)]; 




Alrr 

^ [((Yrl-(ml*ul))*Nvdotl)+(Nrl*(ml-Yvdotl))’ 

]/ 

[{(ml- 

Yvdot1) 

*(Izl-Nrdotl)) 

-(Nvdotl*Yrdotl)]; 




Blr 

= [(T*Nvdotl)-(T*xpl*(ml-Yvdotl))’ 

]/ 

[((ml- 

Yvdotl) 

*(Izl-Nrdotl)) 

-(Nvdotl*Yrdotl)]; 




Clr 

= [(Ydel*Nvdotl)+(Ndel*(ml-Yvdotl))’ 

]/ 

[((ml- 

Yvdotl) 

*(Izl-Nrdotl)) 

g. 

-(Nvdotl*Yrdotl)]; 




0 

% Find 

g. 

control gains 




0 

c 

= zeros(4); 




C(l,3) 

= 1; 




C(2,2) 

= Alvv; 




C(2,3) 

= Alvr; 




C(3,2) 

= Alrv; 




C(3,3) 

= Alrr; 




C(4,l) 

= 1; 




C(4,2) 

= 1; 




D 

= zeros (4,1); 




D(2,l) 

= Civ; 




D(3,l) 

= Clr; 
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poles = [bpole bpole-0.05 bpole-0.10 bpole-0.15]; 
k = place(C,D, poles) ; 

Kp s 1 = k (1, 1) ; 

Kv = k (1,2 ) ; 

Kr = k(l,3) ; 

Ky = k(l,4); 

g, 

0 

% A matrix 

g, 

0 

A = zeros(8); 

A (1,1) = A2vv; 

A (1,2) = A2vr; 

A(1,5) = B2v/L; 

A(1,6) = -B2v/L; 

A(l,7) = B2v+((B2v*xp2)/L); 

A(l,8) = (B2v*xpl)/L; 

A (2,1) = A2rv; 

A (2,2) = A2rr; 

A (2,5) = B2r/L; 

A (2,6) = -B2r/L; 

A(2,7) = B2r+((B2r*xp2)/L) ; 

A(2,8) = (B2r*xpl)/L; 

A(3,3) = Alvv-(Clv*Kv); 

A(3,4) = Alvr-(Clv*Kr) ; 

A(3,5) = Blv/L; 

A(3,6) = -Blv/L-(Clv*Ky); 

A(3,7) = (Blv*xp2)/L; 

A(3,8) = Blv+((Blv*xpl)/L)-(Clv*Kpsi); 

A(4,3) = Alrv-(Clr*Kv); 

A(4,4) = Alrr-(Clr*Kr) ; 

A(4,5) = Blr/L; 

A(4,6) = -Blr/L-(Clr*Ky); 

A(4,7) = (Blr*xp2)/L; 

A(4,8) = Blr+((Blr*xpl)/L)-(Clr*Kpsi); 

A(5,l) = 1; 

A(5,7) = 1; 

A(6,3) = 1; 

A(6,8) = 1; 

A(7,2) = 1; 

A(8,4) = 1; 

"6 

% Find eigenvalues 

g, 

0 

B=eig(A); 

F=real(B); 

H_new=max(F); 

g, 

0 

% Detect if H changes its sign 

g, 

0 

if indexTC==l 

H_old=H_new; 

else 

PROD=H_new*H_old; 
if PROD > 0 

g, 

0 

% H does not change sign - keep going 

33 



H_old=H_new; 

else 

g, 

0 

% H changed its sign - find critical TC by linear 

interpolation 

g, 

0 

index=index+l; 

TC_crit(index)=-(H_old*TC_v(indexTC)- 
H_new*TC_v(indexTC-1))/(H_new-H_old); 

H_old=H_new; 

end 

end 

end 

end 

"0 

% Plot (L, crit TC) curve 
"0 

%plot(L_v,TC_crit),xlabel('L'),ylabel('{T_C}_{\rm 
crit}'),title ('Constant T'),grid 
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APPENDIX B. MATLAB CODE TO FIND THE NONLINEAR 
DYNAMICS COEFFICIENT 


% Nonlinear Dynamics coefficient graph - K vs. L - constant T 


% Parameters: 
% T 
% L 
% TC 
% xpl 
amidships) 

% xp2 = 

amidships) 


Nondimensional tension 

Towline length 

Control law time constant 

Attachment point (leading ship - positive forward of 
Attachment point (trailing ship - positive aft of 


% Constants 


ul 

= 

1; 

u2 

= 

1; 

ml 

= 

0.018078; 

m2 

= 

0.018; 

Izl 

= 

0.0007; 

Iz2 

= 

0.00069412; 

Yvl 

= 

-0.07893; 

Yv2 

= 

-0.1183; 

Yrl 

= 

-0.004044; 

Yr2 

= 

-0.0042; 

Nvl 

= 

-0.016428; 

Nv2 

= 

-0.0187; 

Nrl 

= 

-0.010332; 

Nr2 

= 

-0.0176; 

Yvdotl 

= 

-0.051328; 

Yvdot2 

= 

-0.0184; 

Yrdotl 

= 

0.005617; 

Yrdot2 

= 

-0.0011; 

Nvdotl 

= 

-0.001945; 

Nvdot2 

= 

-0.0008489; 

Nrdotl 

= 

-0.00564; 

Nrdot2 

= 

-0.0090; 

Ydel 

= 

0.0103; 

Ndel 

= 

-0.0051; 

indexL 

= 

0; 

index 

= 

0; 

xpl 

= 

0.5; 

xp2 

= 

0.5; 

deltasat^ 

o 

II 

'0 

% Enter ^ 

constant T 

0 

T 

= 

input('T = 


% Start loop on length 


for iL =1.53:0.01:2.0; 
indexL = indexL+1; 
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indexTC = 

0; 

L 

iL 

L_v(indexL)= 

L; 

A3 

1/L; 

B3 

xpl/L; 

C3 

xp2/L; 

D3 

-1/L; 


% Loop on TC 

for iTC =0.1:0.005 : 2.0; 
indexTC=indexTC+l; 

TC =iTC; 

bpole = -1/TC; 

TC_v(indexTC)=TC; 

% Setup the matrix coefficients 

g, 

0 

A2vv = [((Iz2-Nrdot2)*Yv2)+(Nv2*Yrdot2)]/[((m2-Yvdot2)*(Iz2-Nrdot2))- 
(Nvdot2*Yrdot2)]; 

A2vr = [((Yr2-(m2*u2))*(Iz2-Nrdot2))+(Nr2*Yrdot2)]/[((m2-Yvdot2)*(Iz2- 
Nrdot2))-(Nvdot2*Yrdot2)]; 

B2v = [-((Iz2-Nrdot2)*T)-(Yrdot2*T*xp2)]/[((m2-Yvdot2)*(Iz2-Nrdot2))- 
(Nvdot2*Yrdot2)] ; 

A2rv = [(Yv2*Nvdot2)+(Nv2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2-Nrdot2))- 
(Nvdot2*Yrdot2)]; 

A2rr = [((Yr2-(m2*u2))*Nvdot2)+(Nr2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2- 
Nrdot2))-(Nvdot2*Yrdot2)]; 

B2r = [-(T*Nvdot2)-(T*xp2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2-Nrdot2))- 
(Nvdot2*Yrdot2)]; 

Alvv = [((Izl-Nrdotl)*Yvl)+(Nvl*Yrdotl)]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)]; 

Alvr = [(Yrl-ml*ul)*(Izl-Nrdotl)+(Nrl*Yrdotl)]/[((ml-Yvdotl)*(Izl- 
Nrdotl) )-(Nvdotl*Yrdotl)]; 

Blv = [((Izl-Nrdotl)*T)-(Yrdotl*T*xpl)]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)]; 

Civ = [((Izl-Nrdotl)*Ydel)+(Ndel*Yrdot1)]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)]; 

Alrv = [(Yvl*Nvdotl)+(Nvl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)]; 

Alrr = [((Yrl-(ml*ul))*Nvdotl)+(Nrl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 
Nrdotl) )-(Nvdotl*Yrdotl)]; 

Blr = [(T*Nvdotl)-(T*xpl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)]; 

Clr = [(Ydel*Nvdotl)+(Ndel*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)]; 

% Find control gains 

g, 

0 


C 

= zeros(4) ; 

C(l,3) 

= 1; 

C(2,2) 

= Alvv; 

C(2,3) 

= Alvr; 

C(3,2) 

= Alrv; 

C(3,3) 

= Alrr; 

C(4,l) 

= 1; 

C(4,2) 

= 1; 

D 

= zeros(4,1); 

D(2,l) 

= Civ; 
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D(3,l) 

= 

Clr; 

poles 


[bpole bpole-0.05 

k 


place(C,D,poles) ; 

kpsi 


k(l,l) ; 

kv 


k(l,2) ; 

kr 


k(l,3) ; 

ky 

= 

k(l,4) ; 

% A matrix 

0 

A 


zeros(8); 

A(l, 1) 


A2vv; 

A(l,2) 


A2vr; 

A(l,5) 


B2v/L; 

A(l, 6) 


-B2v/L; 

A(l,7) 


B2v+((B2v*xp2)/L) 

A(l,8) 


(B2v*xpl)/L; 

A(2, 1) 


A2rv; 

A(2,2) 


A2rr; 

A(2,5) 


B2r/L; 

A(2, 6) 


-B2r/L; 

A(2,7) 


B2r+((B2r*xp2)/L) 

A(2,8) 


(B2r*xpl)/L; 

A(3,3) 


Alvv-(Clv*kv); 

A(3,4) 


Alvr-(Clv*kr); 

A(3,5) 


Blv/L; 

A(3, 6) 


-Blv/L-(Clv*ky); 

A(3,7) 


(Blv*xp2)/L; 

A(3,8) 


Blv+((Blv*xpl)/L) 

A(4,3) 


Alrv-(Clr*kv); 

A(4,4) 


Alrr- (Clr*kr); 

A(4,5) 


Blr/L; 

A(4, 6) 


-Blr/L-(Clr*ky); 

A(4,7) 


(Blr*xp2)/L; 

A(4,8) 


Blr+((Blr*xpl)/L) 

A(5, 1) 


1; 

A(5,7) 


1; 

A(6,3) 


1; 

A(6,8) 


1; 

A(7,2) 


1; 

A(8,4) 


1; 


% Find eigenvalues 
"0 

B=eig(A); 

F=real(B); 

H_new=max(F); 

"6 

% Detect if H changes its sign 


if indexTC==l 

H_old=H_new; 

else 

PR0D=H_new*H_old; 
if PROD > 0 
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% H does not change sign - keep going 
"0 

H_old=H_new; 

else 

"6 

% H changed its sign - find critical TC by linear interpolation 
"0 

index=index+l; 

TC_crit(index)=-(H_old*TC_v(indexTC)-H_new*TC_v(indexTC-1))/(H_new- 
H_old) ; 

H_old=H_new; 

"0 

bpole_crit = -l/TC_crit(index); 

poles = [bpole_crit bpole_crit-0.05 bpole_crit-0.10 bpole_crit-0.15]; 

k = place(C,D,poles); 

kpsi = k (1,1) ; 

kv = k (1,2) ; 

kr = k(l,3) ; 

ky = k(l,4) ; 

A_crit = zeros (8); 

A_crit(l,l) = A2vv; 

A_crit(l,2) = A2vr; 

A_crit(l,5) = B2v/L; 

A_crit(l,6) = -B2v/L; 

A_crit(l,7) = B2v+((B2v*xp2)/L); 

A_crit(l,8) = (B2v*xpl)/L; 

A_crit(2,l) = A2rv; 

A_crit(2,2) = A2rr; 

A_crit(2,5) = B2r/L; 

A_crit(2,6) = -B2r/L; 

A_crit(2,7) = B2r+( (B2r*xp2)/L); 

A_crit(2,8) = (B2r*xpl)/L; 

A_crit(3,3) = Alvv-(Clv*kv); 

A_crit(3,4) = Alvr-(Clv*kr); 

A_crit(3,5) = Blv/L; 

A_crit(3,6) = -Blv/L-(Clv*ky); 

A_crit(3,7) = (Blv*xp2)/L; 

A_crit(3,8) = Blv+((Blv*xpl)/L)-(Clv*kpsi); 

A_crit(4,3) = Alrv-(Clr*kv); 

A_crit(4,4) = Alrr-(Clr*kr); 

A_crit(4,5) = Blr/L; 

A_crit(4,6) = -Blr/L-(Clr*ky); 

A_crit(4,7) = (Blr*xp2)/L; 

A_crit(4,8) = Blr+((Blr*xpl)/L)-(Clr*kpsi); 

A_crit(5,1) = 1; 

A_crit(5,7) = 1; 

A_crit(6,3) = 1; 

A_crit(6,8) = 1; 

A_crit(7,2) = 1; 

A_crit(8,4) = 1; 

"6 

%Repeated coloumns are 4-5 and 6-7 of the Devcs 
J=eig(A_crit); 

[D_evcs,V_evls]=eig(A_crit); 
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T_matrixl(:,1)=D_evcs {:,!)} 

T_matrixl(:,2)=D_evcs(:,2); 

T_matrixl(:,3)=D_evcs(:,3); 

T_matrixl(:,4)=real(D_evcs(:,4)); 

T_matrixl(:,5)=imag(D_evcs(: , 4)); 

T_matrixl(:,6)=real(D_evcs(:,6)); 

T_matrixl(:,7)=imag(D_evcs(:,6)); 

T_matrixl(:,8)=D_evcs(:,8); 

"0 

T_matrix(:,1)=T_matrixl(: , 6); 

T_matrix(:,2)=T_matrixl(:,7); 

T_matrix(:,3)=T_matrixl(:,1); 

T_matrix(:,4)=T_matrixl(:,2); 

T_matrix(:,5)=T_matrixl(: , 3); 

T_matrix(:,6)=T_matrixl(: , 4); 

T_matrix (■.,!) =T_matrixl ( : , 5) ; 

T_matrix(:,8)=T_matrixl(: , 8); 

"0 

mll=T_matrix(1,1);ml2=T_matrix(1,2); 
m21=T_matrix(2,1);m22=T_matrix(2,2); 
m31=T_matrix(3,1);m32=T_matrix(3,2) ; 
m4l=T_matrix(4,1);m42=T_matrix(4,2); 
m51=T_matrix(5,1);m52=T_matrix(5,2); 
m8l=T_matrix(8,1);m82=T_matrix(8,2); 
m61=T_matrix(6,1);m62=T_matrix(6,2); 

"6 

T_inverse=inv(T_matrix); 

N1l=T_inverse(1,1) ; 

N12=T_inverse(1,2) ; 

N13=T_inverse (1,3); 

N14=T_inverse (1,4); 

N15=T_inverse (1,5); 

N1l=T_inverse (1,1); 

N2l=T_inverse (2,1); 

N22=T_inverse (2,2); 

N23=T_inverse (2,3); 

N24=T_inverse (2,4); 

N25=T_inverse (2,5); 

"0 

M=inv(T_matrix)*A_crit*T_matrix; 

Lll_l = l/6*B2v/L"3*m41^3-l/6*B2v/L^3*m81"3-l/(6*L"3)*B2v*m51"3*xp2"3; 
Lll_2 = -1/(2*L^2)*B2v*m51"3*xp2"2-l/6*B2v*m51"3-2/(3*L)*B2v*m51"3*xp2; 
Lll_3 = -1/(6*L)*B2v*mll^3*xpl- 

/(6*L"3)*B2v*mll"3*xpl"3+l/(2*L^3)*B2v*m41*m81"2; 

Lll_4 =-1/(2*L^3)*B2v*m41^2*m81-l/L"2*B2v*m81*m51^2*xp2- 
1/(2*L)*B2v*m81*m51^2; 

Lll_5= 1/(2*L"3)*B2v*m81*m51"2*xp2"2+l/L^2*B2v*m41*m51"2*xp2 
+1/(2*L)*B2v*m41*m51^2; 

Lll_6 =+1/(2*L^3)*B2v*m41*m51"2*xp2^2-l/(2*L"3)*B2v*xpl^2*m81*mll"2; 
Lll_7 = -1/(2*L"3)*B2v*xpl^2*m81*mll^2-l/(2*L^3)*B2v*xpl*m81"2*mll; 
Lll_8 = -1/(2*L"2)*B2v*m41^2*m51-l/(2*L"3)*B2v*m41^2*m51*xp2; 

Lll_9 =-1/(2*L^3)*B2v*xpl*m41"2*mll+l/(2*L^3)*B2v*xpl"2*m41*mll"2; 
Lll_10 =+l/L^3*B2v*xp2*m41*m81*m51+l/L^2*B2v*m41*m81*m51 
+l/L^3*B2v*m41*m81*mll*xpl; 
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=-l/L"3*B2v*m81*mll*m51*xp2*xpl-l/L^2*B2v*m81*mll*m51*xpl; 

Lll_12 =+l/L^3*B2v*m41*mll*m51*xp2*xpl+l/L^2*B2v*m41*mll*m51*xpl; 

Lll_13 =-1/(2*L^3)*B2v*mll*m51"2*xp2^2*xpl-l/(2*L)*B2v*mll*m51"2*xpl; 
Lll_14 =-l/L^2*B2v*mll*m51^2*xp2*xpl-l/(2*L^3)*B2v*mll"2*m51*xp2*xpl"2; 
Lll_15 =-1/(2*L^2)*B2v*mll"2*m51*xpl"2; 

Lll = L11_1+L11_2+L11_3+L11_4+L11_5+L11_6+L11_7+L11_8+L11_9+L11_10... 
+L11_11+L11_12+L11_13+L11_14+L11_15; 

L12_l = l/6*B2v/L"3*m42^3-l/6*B2v/L^3*m82"3-l/(6*L^3)*B2v*m52"3*xp2"3; 

L12_2 = -1/(2*L^2)*B2v*m52"3*xp2"2-l/6*B2v*m52"3-2/(3*L)*B2v*m52"3*xp2; 
L12_3 = -1/(6*L)*B2v*ml2^3*xpl- 

1/(6*L^3)*B2v*ml2"3*xpl^3+l/(2*L^3)*B2v*m42*m82^2; 

L12_4 = -1/(2*L^3)*B2v*m42"2*m82-l/L^2*B2v*m82*m52"2*xp2- 
1/(2*L)*B2v*m82*m52"2; 

L12_5 = -1/(2*L^3)*B2v*m82*m52"2*xp2"2+l/L^2*B2v*m42*m52"2*xp2; 

L12_6 = +1/(2*L)*B2v*m42*m52^2+l/(2*L^3)*B2v*m42*m52"2*xp2"2 

L12_7 = -1/(2*L^3)*B2v*xpl"2*m82*ml2^2-l/(2*L"3)*B2v*xpl^2*m82*ml2"2; 

L12_8 = -1/(2*L^3)*B2v*xpl*m82"2*ml2-l/(2*L"2)*B2v*m42"2*m52; 

L12_9 = -1/(2*L"3)*B2v*m42^2*m52*xp2-l/(2*L"3)*B2v*xpl*m42^2*ml2; 

L12_10 = +1/(2*L"3)*B2v*xpl^2*m42*ml2"2+l/L^3*B2v*xp2*m42*m82*m52; 

L12_ll = +l/L^2*B2v*m42*m82*m52+l/L^3*B2v*m42*m82*ml2*xpl; 

L12_12 = -l/L^3*B2v*m82*ml2*m52*xp2*xpl-l/L^2*B2v*m82*ml2*m52*xpl; 

L12_13 = +l/L"3*B2v*m42*ml2*m52*xp2*xpl+l/L^2*B2v*m42*ml2*m52*xpl; 

L12_14 = -1/(2*L^3)*B2v*ml2*m52"2*xp2^2*xpl-l/(2*L)*B2v*ml2*m52"2*xpl; 
L12_15 = -l/L"2*B2v*ml2*m52"2*xp2*xpl-l/(2*L"3)*B2v*ml2^2*m52*xp2*xpl"2; 
L12_16 = -1/(2*L"2)*B2v*ml2^2*m52*xpl"2; 

L12 = L12_l+L12_2+L12_3+L12_4+L12_5+L12_6+L12_7+L12_8+L12_9+L12_10... 

+L12_11+L12_12+L12_13+L12_14+L12_15+L12_16; 

L13_l=l/2*B2v/L^3*m41"2*m42-l/2*B2v/L"3*m81"2*m82 
-1/(2*L^3)*B2v*m51^2*m52*xp2^3; 

L13_2=-3/(2*L"3)*B2v*m51"2*m52*xp2"2-l/2*B2v*m51"2*m52- 
2/L*B2v*m51^2*m52*xp2; 

L13_3 = -1/(2*L)*B2v*mll"2*ml2*xpl-l/(2*L"3)*B2v*mll"2*ml2*xpl"3; 
L13_4=+l/L^3*B2v*m41*m81*m82+l/(2*L^3)*B2v*m42*m81"2- 
1/(2*L^3)*B2v*m41^2*m82; 

L13_5=-l/L^3*B2v*m41*m42*m81-2/(L"2)*B2v*m81*m51*m52*xp2- 
l/L*B2v*m81*m51*m52; 

L13_6 = -1/(L"3)*B2v*m81*m51*m52*xp2"2-l/(L^2)*B2v*m82*m51"2*xp2; 

L13_7 = -1/(2*L)*B2v*m82*m51"2-l/(2*L"3)*B2v*m82*m51"2*xp2"2; 

L13_8=2/(L^2)*B2v*m41*m51*m52*xp2+l/L*B2v*m41*m51*m52+l/L"3*B2v*m41*m51* 
m52*xp2^2; 

L13_9=l/L^2*B2v*m42*m51"2*xp2+l/(2*L)*B2v*m42*m51^2+l/(2*L^3)*B2v*m42*m5 
l"2*xp2^2; 

L13_10 = -l/L"3*B2v*xpl^2*m81*mll*ml2-l/(2*L^3)*B2v*xpl^2*m82*mll"2; 
L13_ll = -l/L"3*B2v*xpl^2*m81*mll*ml2-l/(2*L"3)*B2v*xpl^2*m82*mll"2; 
L13_12 = -1/(2*L"3)*B2v*xpl*m81^2*ml2-l/L"3*B2v*xpl*m81*m82*mll; 
L13_13=-l/(2*L"2)*B2v*m41^2*m52-l/(2*L^2)*B2v*m41"2*m52*xp2- 
l/L^2*B2v*m41*m42*m51; 

L13_14 = -l/L"3*B2v*m41*m42*m51*xp2-l/(2*L"3)*B2v*xpl*m41"2*ml2; 

L13_15 = -l/L^3*B2v*xpl*m41*m42*mll+l/L^3*B2v*xpl"2*m41*mll*ml2; 

L13_16 = +1/(2*L"3)*B2v*xpl^2*m42*mll"2+l/L^2*B2v*xp2*m41*m81*m52; 
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L13_17=l/L^2*B2v*m41*m81*m52+l/L^3*B2v*xp2*m41*m82*m51+l/L^2*B2v*m41*m82 

*m51; 

L13_18=+l/L^3*B2v*xp2*m42*m81*m51+l/L^2*B2v*m42*m81*m51+l/L"3*B2v*m41*m8 
l*ml2*xpl; 

L13_19 =+l/L"3*B2v*m41*m82*mll*xpl+l/L"3*B2v*m42*m81*mll*xpl; 

L13_20 = -l/L^3*B2v*m81*mll*m52*xp2*xpl-l/L"2*B2v*m81*mll*m52*xpl; 

L13_21 = -l/L"3*B2v*m81*ml2*m51*xp2*xpl-l/L^2*B2v*m81*ml2*m51*xpl; 

L13_22 = -l/L"3*B2v*m82*mll*m51*xp2*xpl-l/L^2*B2v*m82*mll*m51*xpl; 

L13_23 = +l/L^3*B2v*m41*mll*m52*xp2*xpl+l/L^2*B2v*m41*mll*m52*xpl; 

L13_24 = +l/L"3*B2v*m41*ml2*m51*xp2*xpl+l/L^2*B2v*m41*ml2*m51*xpl; 

L13_25 = +l/L"3*B2v*m42*mll*m51*xp2*xpl+l/L^2*B2v*m42*mll*m51*xpl; 

L13_26 = -l/L"3*B2v*mll*m51*m52*xp2^2*xpl-l/L*B2v*mll*m51*m52*xpl; 
L13_27=-2/L^2*B2v*mll*m51*m52*xp2*xpl-l/(2*L^3)*B2v*ml2*m51"2*xp2"2*xpl; 
L13_28 = -1/(2*L)*B2v*ml2*m51^2*xpl-l/L^2*B2v*ml2*m51^2*xp2*xpl; 
L13_29=-l/(2*L"'3)*B2v*mll"2*m52*xp2*xpl"2-l/(2*L"2)*B2v*mll"2*m52*xpl"2; 
L13_30 = -l/L^3*B2v*mll*ml2*m51*xp2*xpl"2-l/L^2*B2v*mll*ml2*m51*xpl"2; 

L13 = L13_l+L13_2+L13_3+L13_4+L13_5+L13_6+L13_7+L13_8+L13_9+L13_10... 

+L13_11+L13_12+L13_13+L13_14+L13_15+L13_16+L13_17+L13_18+L13_190 
+L13_20+L13_21+L13_22+L13_23+L13_24+L13_25+L13_26+L13_27+L13_28O 
+L13_29+L13_30; 

L14_l=l/2*B2v/L^3*m41*m42"2-l/2*B2v/L^3*m81*m82^2- 
1/(2*L^3)*B2v*m51*m52^2*xp2^3; 

L14_2=-3/(2*L"2)*B2v*m51*m52"2*xp2"2-l/2*B2v*m51*m52"2- 
2/L*B2v*m51*m52^2*xp2; 

L14_3=-l/(2*L)*B2v*mll*ml2"2*xpl-l/(2*L^3)*B2v*mll*ml2"2*xpl"3; 

L14_4=+l/(2*L"3)*B2v*m41*m82"2+l/L"3*B2v*m42*m81*m82- 
l/L^3*B2v*m41*m42*m82; 

L14_5=-l/(2*L"3)*B2v*m42"2*m81-l/(L^2)*B2v*m81*m52"2*xp2- 
1/(2*L)*B2v*m81*m52"2; 

L14_6=-l/(2*L"3)*B2v*m81*m52"2*xp2"2-2/(L^2)*B2v*m82*m51*m52*xp2; 
L14_7=-l/L*B2v*m82*m51*m52- 

l/L^3*B2v*m82*m51*m52*xp2"2+l/L^2*B2v*m41*m52^2*xp2; 

L14_8=+l/(2*L)*B2v*m41*m52"2*+l/(2*L"3)*B2v*m41*m52^2*xp2"2; 

L14_9=+2/(L^2)*B2v*m42*m51*m52*xp2+l/L*B2v*m42*m51*m52+l/L^3*B2v*m42*m51 
*m52*xp2^2; 

L14_10=-l/(2*L"3)*B2v*xpl"2*m81*ml2^2-l/L^3*B2v*xpl^2*m82*mll*ml2; 
L14_ll=-1/(2*L^3)*B2v*xpl^2*m81*ml2^2-l/L^3*B2v*xpl^2*m82*mll*ml2; 
L14_12=-l/L"3*B2v*xpl*m81*m82*ml2-l/(2*L"3)*B2v*xpl*m82^2*mll; 
L14_13=-l/L"2*B2v*m41*m42*m52-l/L^3*B2v*m41*m42*m52*xp2- 
1/(2*L"2)*B2v*m42"2*m51; 

L14_14=-l/(2*L"3)*B2v*m42"2*m51*xp2-l/L^3*B2v*xpl*m41*m42*ml2; 

L14_15=-l/(2*L^3)*B2v*xpl*m42^2*mll+l/(2*L^3)*B2v*xpl^2*m41*ml2"2; 
L14_16=+l/L^3*B2v*xpl"2*m42*mll*ml2+l/L^3*B2v*xp2*m41*m82*m52+l/L"2*B2v* 
m41*m82*m52; 

L14_17=+l/L^3*B2v*xp2*m42*m81*m52+l/L"2*B2v*m42*m81*m52+l/L"3*B2v*xp2*m4 

2*m82*m51; 

L14_18=+l/L^2*B2v*m42*m82*m51+l/L^3*B2v*m41*m82*ml2*xpl+l/L"3*B2v*m42*m8 
l*ml2*xpl; 

L14_19=+l/L^3*B2v*m42*m82*mll*xpl-l/L^3*B2v*m81*ml2*m52*xp2*xpl; 

L14_20=-l/L"2*B2v*m81*ml2*m52*xpl-l/L^3*B2v*m82*mll*m52*xp2*xpl; 

L14_21=-l/L^2*B2v*m82*mll*m52*xpl-l/L"3*B2v*m82*ml2*m51*xp2*xpl; 

L14_22=-l/L^2*B2v*m82*ml2*m51*xpl+l/L^3*B2v*m41*ml2*m52*xp2*xpl; 

L14_23=+l/L^2*B2v*m41*ml2*m52*xpl+l/L"3*B2v*m42*mll*m52*xp2*xpl; 
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L14_24=+l/L^2*B2v*m42*mll*m52*xpl+l/L^3*B2v*m42*ml2*m51*xp2*xpl; 
L14_25=+l/L^2*B2v*m42*ml2*m51*xpl-l/(2*L"3)*B2v*mll*m52^2*xp2"2*xpl; 
L14_26=-l/(2*L)*B2v*mll*m52^2*xpl-l/L^2*B2v*mll*m52^2*xp2*xpl; 

L14_27=-l/(L^3)*B2v*ml2*m51*m52*xp2^2*xpl-l/L*B2v*ml2*m51*m52*xpl; 
L14_28=-2/L"2*B2v*ml2*m51*m52*xp2*xpl-l/(L"3)*B2v*mll*ml2*m52*xp2*xpl"2; 
L14_29=-l/(L^2)*B2v*mll*ml2*m52*xpl"2-l/(2*L"3)*B2v*ml2"2*m51*xp2*xpl"2; 
L14_30= -l/L"2*B2v*ml2"2*m51*xpl"2; 

L14=L14_1+L14_2+L14_3+L14_4+L14_5+L14_6+L14_7+L14_8+L14_9+L14_10... 

+L14_11+L14_12+L14_13+L14_14+L14_15+L14_16+L14_17+L14_18+L14_19+L14_20.. 

+L14_21+L14_22+L14_23+L14_24+L14_25+L14_26+L14_27+L14_28+L14_29+L14_30; 

L21_l=l/(6*L^3)*B2r*m41^3-l/(6*L"3)*B2r*m81^3-l/(6*L"3)*B2r*m51"3*xp2"3; 
L21_2=-l/(2*L^2)*B2r*m51"3*xp2^2-l/6*B2r*m51^3-2/(3*L)*B2r*m51"3*xp2; 
L21_3=-l/(6*L)*B2r*mll"3*xpl- 

1/(6*L^3)*B2r*mll^3*xpl^3+l/(2*L^3)*B2r*m41*m81^2; 

L21_4=-l/(2*L^3)*B2r*m41^2*m81-l/(L^2)*B2r*m81*m51^2*xp2- 
1/(2*L)*B2r*m81*m51^2; 

L21_5=-l/(2*L^3)*B2r*m81*m51"2*xp2^2+l/(L^2)*B2r*m41*m51"2*xp2; 

L21_6=+l/(2*L)*B2r*m41*m51"2+l/(2*L^3)*B2r*m41*m51"2*xp2"2; 

L21_7=-l/(2*L^3)*B2r*xpl^2*m81*mll^2-l/(2*L^3)*B2r*xpl^2*m81*mll"2; 
L21_8=-l/(2*L^3)*B2r*xpl*m81^2*mll-l/(2*L"2)*B2r*m41^2*m51; 

L21_9=-l/(2*L"3)*B2r*m41"2*m51*xp2-l/(2*L"3)*B2r*xpl*m41"2*mll; 
L21_10=+l/(2*L^3)*B2r*xpl^2*m41*mll"2+l/(L^3)*B2r*xp2*m41*m81*m51; 
L21_ll=+1/(L"2)*B2r*m41*m81*m51+l/(L"3)*B2r*m41*m81*mll*xpl; 

L21_12=-l/(L^3)*B2r*m81*mll*m51*xp2*xpl-l/(L"2)*B2r*m81*mll*m51*xpl; 
L21_13=+l/(L^3)*B2r*m41*mll*m51*xp2*xpl+l/(L"2)*B2r*m41*mll*m51*xpl; 
L21_14=-l/(2*L"3)*B2r*mll*m51"2*xp2^2*xpl-l/(2*L)*B2r*mll*m51"2*xpl; 
L21_15=-l/(L^2)*B2r*mll*m51"2*xp2*xpl-l/(2*L"3)*B2r*mll^2*m51*xp2*xpl"2; 
L21_16=-l/(2*L^2)*B2r*mll"2*m51*xpl^2; 

L21= L21_l+L21_2+L21_3+L21_4+L21_5+L21_6+L21_7+L21_8+L21_9+L21_10... 

+L21_11+L21_12+L21_13+L21_14+L21_15+L21_16; 

L22_l=l/(6*L"3)*B2r*m42^3-l/(6*L^3)*B2r*m82^3-l/(6*L"3)*B2r*m52"3*xp2"3; 
L22_2=-l/(2*L"2)*B2r*m52"3*xp2"2-l/6*B2r*m52"3-2/(3*L)*B2r*m52"3*xp2; 
L22_3=-l/(6*L)*B2r*ml2"3*xpl- 

1/(6*L^3)*B2r*ml2"3*xpl^3+l/(2*L^3)*B2r*m42*m82^2; 

L22_4=-l/(2*L^3)*B2r*m42^2*m82-l/(L"2)*B2r*m82*m52^2*xp2- 
1/(2*L)*B2r*m82*m52^2; 

L22_5=-l/(2*L"3)*B2r*m82*m52"2*xp2"2+l/(L^2)*B2r*m42*m52"2*xp2; 

L22_6=+l/(2*L)*B2r*m42*m52"2+l/(2*L^3)*B2r*m42*m52^2*xp2^2; 

L22_7=-l/(2*L"3)*B2r*xpl"2*m82*ml2"2-l/(2*L^3)*B2r*xpl"2*m82*ml2"2; 
L22_8=-l/(2*L^3)*B2r*xpl*m82^2*ml2-l/(2*L^2)*B2r*m42^2*m52; 

L22_9=-l/(2*L"3)*B2r*m42"2*m52*xp2-l/(2*L"3)*B2r*xpl*m42"2*ml2; 
L22_10=+l/(2*L"3)*B2r*xpl"2*m42*ml2^2+l/(L"3)*B2r*xp2*m42*m82*m52; 
L22_ll=+1/(L"2)*B2r*m42*m82*m52+l/(L^3)*B2r*m42*m82*ml2*xpl; 

L22_12=-l/(L^3)*B2r*m82*ml2*m52*xp2*xpl-l/(L"2)*B2r*m82*ml2*m52*xpl; 
L22_13= +1/(L^3)*B2r*m42*ml2*m52*xp2*xpl+l/(L^2)*B2r*m42*ml2*m52*xpl; 
L22_14=-l/(2*L^3)*B2r*ml2*m52^2*xp2"2*xpl-l/(2*L)*B2r*ml2*m52"2*xpl; 
L22_15 =-1/(L"2)*B2r*ml2*m52"2*xp2*xpl- 
1/(2*L^3)*B2r*ml2"2*m52*xp2*xpl^2; 

L22_16 = -1/(2*L"2)*B2r*ml2^2*m52*xpl"2; 
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L22 


= L22_l+L22_2+L22_3+L22_4+L22_5+L22_6+L22_7+L22_8+L22_9+L22_10... 
+L22_11+L22_12+L22_13+L22_14+L22_15+L22_16; 

L23_l=l/(2*L^3)*B2r*m41"2*m42-l/(2*L^3)*B2r*m81"2*m82- 
1/(2*L"3)*B2r*m51"2*m52*xp2^3; 

L23_2=-3/(2*L"2)*B2r*m51^2*m52*xp2^2-l/2*B2r*m51^2*m52- 
2/L*B2r*m51^2*m52*xp2; 

L23_3=-l/(2*L)*B2r*mll"2*ml2*xpl-l/(2*L^3)*B2r*mll"2*ml2*xpl"3; 

L23_4=+l/(L"3)*B2r*m41*m81*m82+l/(2*L^3)*B2r*m42*m81"2- 
1/(2*L"3)*B2r*m41^2*m82 ; 

L23_5=-l/(L^3)*B2r*m41*m42*m81-2/(L^2)*B2r*m81*m51*m52*xp2- 
l/L*B2r*m81*m51*m52; 

L23_6=-l/(L"3)*B2r*m81*m51*m52*xp2"2-l/(2*L^2)*B2r*m82*m51"2*xp2; 
L23_7=-l/(2*L)*B2r*m82*m51^2-l/(2*L"3)*B2r*m82*m51^2*xp2"2; 

L23_8=+2/(L^2)*B2r*m41*m51*m52*xp2+l/L*B2r*m41*m51*m52; 

L23_9=+l/(L^3)*B2r*m41*m51*m52*xp2^2+l/(L^2)*B2r*m42*m51^2*xp2; 
L23_10=+l/(2*L)*B2r*m42*m51^2+l/(2*L^3)*B2r*m42*m51^2*xp2^2; 

L23_ll=-1/(L^3)*B2r*xpl^2*m81*mll*ml2-l/(2*L"3)*B2r*xpl^2*m82*mll"2; 
L23_12=-l/(L^3)*B2r*xpl"2*m81*mll*ml2-l/(2*L^3)*B2r*xpl"2*m82*mll"2; 
L23_13=-l/(2*L^3)*B2r*xpl*m81"2*ml2-l/(L"3)*B2r*xpl*m81*m82*mll; 
L23_14=-l/(2*L^2)*B2r*m41"2*m52-l/(2*L^3)*B2r*m41^2*m52*xp2; 

L23_15=-l/(L^2)*B2r*m41*m42*m51-l/(L^3)*B2r*m41*m42*m51*xp2; 

L23_16=-l/(2*L"3)*B2r*xpl*m41"2*ml2-l/(L"3)*B2r*xpl*m41*m42*mll; 
L23_17=+l/(L^3)*B2r*xpl"2*m41*mll*ml2+l/(2*L^3)*B2r*xpl"2*m42*mll"2; 
L23_18=+l/(L"3)*B2r*xp2*m41*m81*m52+l/(L"2)*B2r*m41*m81*m52; 

L23_19=+l/(L"3)*B2r*xp2*m41*m82*m51+l/(L"2)*B2r*m41*m82*m51; 

L23_20=+l/(L^3)*B2r*xp2*m42*m81*m51+l/(L^2)*B2r*m42*m81*m51; 

L23_21=+l/(L^3)*B2r*m41*m81*ml2*xpl+l/(L"3)*B2r*m41*m82*mll*xpl; 
L23_22=+l/(L^3)*B2r*m42*m81*mll*xpl-l/(L^3)*B2r*m81*mll*m52*xp2*xpl; 
L23_23=-l/(L^2)*B2r*m81*mll*m52*xpl-l/(L"3)*B2r*m81*ml2*m51*xp2*xpl; 
L23_24=-l/(L"2)*B2r*m81*ml2*m51*xpl-l/(L"3)*B2r*m82*mll*m51*xp2*xpl; 
L23_25=-l/(L^2)*B2r*m82*mll*m51*xpl+l/(L^3)*B2r*m41*mll*m52*xp2*xpl; 
L23_26=+l/(L"2)*B2r*m41*mll*m52*xpl+l/(L"3)*B2r*m41*ml2*m51*xp2*xpl; 
L23_27=+l/(L^2)*B2r*m41*ml2*m51*xpl+l/(L^3)*B2r*m42*mll*m51*xp2*xpl; 
L23_28=+l/(L"2)*B2r*m42*mll*m51*xpl-l/(L^3)*B2r*mll*m51*m52*xp2"2*xpl; 
L23_29=-l/L*B2r*mll*m51*m52*xpl-2/(L"2)*B2r*mll*m51*m52*xp2*xpl; 
L23_30=-l/(2*L^3)*B2r*ml2*m51^2*xp2^2*xpl-l/(2*L)*B2r*ml2*m51"2*xpl; 
L23_31=-l/(L"2)*B2r*ml2*m51^2*xp2*xpl-l/(2*L"3)*B2r*mll^2*m52*xp2*xpl"2; 
L23_32=-l/(2*L^2)*B2r*mll^2*m52*xpl"2-l/(L^3)*B2r*mll*ml2*m51*xp2*xpl"2; 
L23_33=-l/(L^2)*B2r*mll*ml2*m51*xpl^2; 


L23=L23_1+L23_2+L23_3+L23_4+L23_5+L23_6+L23_7+L23_8+L23_9+L23_10... 

+L23_11+L23_12+L23_13+L23_14+L23_15+L23_16+L23_17+L23_18+L23_19+L23_20.. 

+L23_21+L23_22+L23_23+L23_24+L23_25+L23_26+L23_27+L23_28+L23_29+L23_30.. 
+L23_31+L23_32+L23_33; 


L24_l=l/(2*L"3)*B2r*m41*m42^2-l/(2*L"3)*B2r*m81*m82^2- 
1/(2*L^3)*B2r*m51*m52^2*xp2^3; 
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L24_2=-3/(2*L"2)*B2r*m51*m52^2*xp2^2-l/2*B2r*m51*m52^2- 
2/L*B2r*m51*m52^2*xp2; 

L24_3 = -1/(2*L)*B2r*mll*ml2^2*xpl-l/(2*L"3)*B2r*mll*ml2^2*xpl"3; 
L24_4=+l/(2*L^3)*B2r*m41*m82^2+l/(L^3)*B2r*m42*m81*m82- 
1/(L"3)*B2r*m41*m42*m82; 

L24_5=-l/(2*L^3)*B2r*m42"2*m81-l/(L"2)*B2r*m81*m52^2*xp2- 
1/(2*L)*B2r*m81*m52^2; 

L24_6=-l/(2*L"3)*B2r*m81*m52"2*xp2^2-l/(2*L^2)*B2r*m82*m51*m52*xp2; 
L24_7=-l/L*B2r*m82*m51*m52- 

1/(L"3)*B2r*m82*m51*m52*xp2^2+l/(L"2)*B2r*m41*m52"2*xp2; 

L24_8=+l/ (2*L) *B2r*m41*m52^2 + l/ (2*L"3) *B2r*m41*m52"'2*xp2^2; 

L24_9=+2/(L^2)*B2r*m42*m51*m52*xp2+l/L*B2r*m42*m51*m52+l/(L"3)*B2r*m42*m 
51*m52*xp2"2; 

L24_10=-l/(2*L^3)*B2r*xpl"2*m81*ml2^2-l/(L^3)*B2r*xpl^2*m82*mll*ml2; 
L24_ll=-1/(2*L"'3)*B2r*xpl"2*m81*ml2^2-l/(L"3)*B2r*xpl"2*m82*mll*ml2; 
L24_12=-l/(L^3)*B2r*xpl*m81*m82*ml2-l/(2*L^3)*B2r*xpl*m82^2*mll; 
L24_13=-l/(L"2)*B2r*m41*m42*m52-l/(L"3)*B2r*m41*m42*m52*xp2- 
1/(2*L^2)*B2r*m42^2*m51; 

L24_14=-l/(2*L^3)*B2r*m42^2*m51*xp2-l/(L^3)*B2r*xpl*m41*m42*ml2; 
L24_15=-l/(2*L^3)*B2r*xpl*m42"2*mll+l/(2*L"3)*B2r*xpl"2*m41*ml2"2; 
L24_16=+l/(L^3)*B2r*xpl"2*m42*mll*ml2+l/(L^3)*B2r*xp2*m41*m82*m52; 
L24_17=+l/(L"2)*B2r*m41*m82*m52+l/(L"3)*B2r*xp2*m42*m81*m52+l/(L"2)*B2r* 
m42*m81*m52; 

L24_18=+l/(L^3)*B2r*xp2*m42*m82*m51+l/(L^2)*B2r*m42*m82*m51; 

L24_19=+l/(L"3)*B2r*m41*m82*ml2*xpl+l/(L"3)*B2r*m42*m81*ml2*xpl; 
L24_20=+l/(L^3)*B2r*m42*m82*mll*xpl-l/(L^3)*B2r*m81*ml2*m52*xp2*xpl; 
L24_21=-l/(L"2)*B2r*m81*ml2*m52*xpl-l/(L"3)*B2r*m82*mll*m52*xp2*xpl; 
L24_22=-l/(L"2)*B2r*m82*mll*m52*xpl-l/(L^3)*B2r*m82*ml2*m51*xp2*xpl; 
L24_23=-l/(L^2)*B2r*m82*ml2*m51*xpl+l/(L^3)*B2r*m41*ml2*m52*xp2*xpl; 
L24_24=+l/(L^2)*B2r*m41*ml2*m52*xpl+l/(L^3)*B2r*m42*mll*m52*xp2*xpl; 
L24_25=+l/(L^2)*B2r*m42*mll*m52*xpl+l/(L^3)*B2r*m42*ml2*m51*xp2*xpl; 
L24_26=+l/(L"2)*B2r*m42*ml2*m51*xpl-l/(2*L"3)*B2r*mll*m52^2*xp2"2*xpl; 
L24_27=-l/(2*L)*B2r*ml1*m52^2*xpl-l/(L"2)*B2r*mll*m52"2*xp2*xpl; 
L24_28=-l/(L^3)*B2r*ml2*m51*m52*xp2"2*xpl-l/L*B2r*ml2*m51*m52*xpl; 
L24_29=-2/(L^2)*B2r*ml2*m51*m52*xp2*xpl- 
1/(L^3)*B2r*mll*ml2*m52*xp2*xpl"2; 

L24_30=-l/(L"2)*B2r*mll*ml2*m52*xpl^2-l/(2*L"3)*B2r*ml2^2*m51*xp2*xpl"2; 
L24_31=-l/(2*L"2)*B2r*ml2"2*m51*xpl^2; 

L24= L24_l+L24_2+L24_3+L24_4+L24_5+L24_6+L24_7+L24_8+L24_9+L24_10... 

+L24_11+L24_12+L24_13+L24_14+L24_15+L24_16+L24_17+L24_18+L24_19+L24_20.. 

+L24_21+L24_22+L24_23+L24_24+L24_25+L24_26+L24_27+L24_28+L24_29+L24_30.. 
+L24_31; 

L31_l=l/(6*L"3)*Blv*m41^3-l/(6*L^3)*Blv*m81^3-l/(6*L)*Blv*m51"3*xp2; 
L31_2=-l/(6*L"3)*Blv*m51"3*xp2^3-l/(2*L^2)*Blv*mll"3*xpl"2- 
l/6*Blv*mll^3; 

L31_3=-2/(3*L)*Blv*mll"3*xpl- 

1/(6*L^3)*Blv*mll^3*xpl"3+l/(2*L^3)*Blv*m41*m81"2; 

L31_4=+l/(2*L"3)*Blv*m41*m81^2-l/(2*L^3)*Blv*xp2^2*m81*m51^2; 

L31_5=+l/(2*L"3)*Blv*xp2"2*m41*m51"2-l/(L^2)*Blv*m81*mll"2*xpl; 

L31_6=-l/(2*L^3)*Blv*m81*mll^2*xpl"2-l/(2*L)*Blv*m81*mll^2; 

L31_7=-l/(2*L"3)*Blv*xp2*m81"2*m51-l/(2*L"2)*Blv*m81"2*mll; 
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L31_8=-l/(2*L^3)*Blv*m81"2*mll*xpl-l/(2*L^3)*Blv*xp2*m41"2*m51; 

L31_9=-l/(2*L"2)*Blv*m41"2*mll- 

1/(2*L^3)*Blv*m41"2*mll*xpl+l/(2*L)*Blv*m41*mll^2; 

L31_10=+l/(L^2)*Blv*m41*mll^2*xpl+l/(2*L^3)*Blv*m41*mll^2*xpl"2; 
L31_ll=+1/(L^3)*Blv*xp2*m41*m81*m51+l/(L"2)*Blv*m41*m81*mll; 

L31_12=+l/(L^3)*Blv*m41*m81*mll*xpl-l/(L^3)*Blv*m81*mll*m51*xp2*xpl; 
L31_13=-l/(L"2)*Blv*m81*mll*m51*xp2+l/(L"3)*Blv*m41*mll*m51*xp2*xpl; 
L31_14=+l/(L^2)*Blv*m41*mll*m51*xp2-l/(2*L^2)*Blv*mll*m51"2*xp2"2; 
L31_15=-l/(2*L^3)*Blv*mll*m51^2*xp2^2*xpl-l/(2*L)*Blv*mll^2*m51*xp2; 
L31_16=-l/(L"2)*Blv*mll^2*m51*xp2*xpl-l/(2*L"3)*Blv*mll^2*m51*xp2*xpl"2; 
L31_17=-l/3*Clv/deltasat^2* k:psi^3*ml1^3-1/3*Clv/deltasat^2*kv^3*m21^3; 
L31_18=-l/3*Clv/deltasat"2*kr"3*m31^3-l/3*Clv/deltasat"2*ky"3*m41"3; 
L31_19=-Clv/deltasat^2*kpsi^2*kv*m21*mll^2- 
Clv/deltasat^2*kpsi^2*kr*m31*ml1^2; 

L31_20=-Clv/deltasat^2*kpsi^2*ky*m41*mll^2- 
Clv/deltasat^2*kpsi*kv^2*m21^2*ml1; 

L31_21=-Clv/deltasat^2*kpsi*kr^2*m31^2*mll- 
Clv/deltasat^2*kpsi*ky^2*m41^2*ml1; 
L31_22=-Clv/deltasat^2*kv^2*kr*m31*m21^2- 
Clv/deltasat^2*kv^2*ky*m4l*m21^2; 
L31_23=-Clv/deltasat^2*kv*kr^2*m21*m31^2- 
Clv/deltasat^2*kv*ky^2*m2l*m41^2; 

L31_24 = -Clv/deltasat^2*kr^2*ky*m41*m31^2- 
Clv/deltasat^2*kr*ky^2*m31*m41^2; 

L31_25 = -2*Clv/deltasat^2*kpsi*kv*kr*m31*m21*ml1; 

L31_26 = -2*Clv/deltasat^2*kpsi*kv*ky*m4l*m2l*ml1; 

L31_27 = -2*Clv/deltasat^2*kpsi*kr*ky*m4l*m31*ml1; 

L31_28 = -2*Clv/deltasat"2*kv*kr*ky*m41*m21*m31; 

L31 = L31_l+L31_2+L31_3+L31_4+L31_5+L31_6+L31_7+L31_8+L31_9+L31_10... 

+L31_11+L31_12+L31_13+L31_14+L31_15+L31_16+L31_17+L31_18+L31_19+L31_20.. 
+L31_21+L31_22+L31_23+L31_24+L31_25+L31_26+L31_27+L31_28; 


L32_l=l/(6*L"3)*Blv*m42^3-l/(6*L"3)*Blv*m82^3-l/(6*L)*Blv*m52"3*xp2; 
L32_2=-l/(6*L"3)*Blv*m52"3*xp2"3-l/(2*L^2)*Blv*ml2"3*xpl"2- 
l/6*Blv*ml2^3; 

L32_3=-2/(3*L)*Blv*ml2"3*xpl- 

1/(6*L^3)*Blv*ml2^3*xpl"3+l/(2*L^3)*Blv*m42*m82"2; 

L32_4=+l/(2*L"3)*Blv*m42*m82"2-l/(2*L"3)*Blv*xp2"2*m82*m52"2; 

L32_5=+l/(2*L"3)*Blv*xp2"2*m42*m52"2-l/(L^2)*Blv*m82*ml2"2*xpl; 

L32_6=-l/(2*L^3)*Blv*m82*ml2^2*xpl^2-l/(2*L)*Blv*m82*ml2^2; 

L32_7=-l/(2*L"3)*Blv*xp2*m82"2*m52-l/(2*L"2)*Blv*m82"2*ml2; 

L32_8=-l/(2*L^3)*Blv*m82"2*ml2*xpl-l/(2*L^3)*Blv*xp2*m42^2*m52; 

L32_9=-l/(2*L"2)*Blv*m42"2*ml2- 

1/(2*L^3)*Blv*m42"2*ml2*xpl+l/(2*L)*Blv*m42*ml2^2; 

L32_10=+l/(L^2)*Blv*m42*ml2^2*xpl+l/(2*L^3)*Blv*m42*ml2"2*xpl"2; 
L32_ll=+1/(L^3)*Blv*xp2*m42*m82*m52+l/(L"2)*Blv*m42*m82*ml2; 

L32_12=+l/(L^3)*Blv*m42*m82*ml2*xpl-l/(L^3)*Blv*m82*ml2*m52*xp2*xpl; 
L32_13=-l/(L^2)*Blv*m82*ml2*m52*xp2+l/(L^3)*Blv*m42*ml2*m52*xp2*xpl; 
L32_14=+l/(L"2)*Blv*m42*ml2*m52*xp2-l/(2*L"2)*Blv*ml2*m52^2*xp2"2; 
L32_15=-l/(2*L"'3)*Blv*ml2*m52^2*xp2^2*xpl-l/(2*L)*Blv*ml2"2*m52*xp2; 
L32_16=-l/(L^2)*Blv*ml2^2*m52*xp2*xpl-l/(2*L"3)*Blv*ml2"2*m52*xp2*xpl"2; 
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L32_17=-l/3*Clv/deltasat^2*k:psi^3*ml2^3-l/3*Clv/deltasat^2*kv^3*m22^3; 
L32_18=-l/3*Clv/deltasat"2*kr"3*m32^3-l/3*Clv/deltasat"2*ky"3*m42"3; 
L32_l9=-Clv/deltasat^2*kpsi^2*kv*m22*ml2^2- 
Clv/deltasat^2*kpsi^2*kr*m32*ml2^2; 

L32_20=-Clv/deltasat^2*kpsi^2*ky*m42*ml2^2- 
Clv/deltasat^2*kpsi*kv^2*m22^2*ml2; 

L32_21=-Clv/deltasat^2*kpsi*kr^2*m32^2*ml2- 
Clv/deltasat^2*kpsi*ky^2*m42^2*ml2; 

L32_22=-Clv/deltasat^2*kv^2*kr*m32*m22^2- 
Clv/deltasat^2*kv^2*ky*m42*m22^2; 

L32_23=-Clv/deltasat^2*kv*kr^2*m22*m32^2- 
Clv/deltasat^2*kv*ky^2*m22*m42^2; 
L32_24=-Clv/deltasat^2*kr"2*ky*m42*m32"2- 
Clv/deltasat^2*kr*ky^2*m32*m42^2; 

L32_25 = -2*Clv/deltasat^2*kpsi*kv*kr*m32*m22*ml2; 

L32_26 = -2*Clv/deltasat^2*kpsi*kv*ky*m42*m22*ml2; 

L32_27 = -2*Clv/deltasat^2*kpsi*kr*ky*m42*m32*ml2; 

L32_28 = -2*Clv/deltasat^2*kv*kr*ky*m42*m22*m32; 


L32 = L32_l+L32_2+L32_3+L32_4+L32_5+L32_6+L32_7+L32_8+L32_9+L32_10... 
+L32_11+L32_12+L32_13+L32_14+L32_15+L32_16+L32_17+L32_18+L32_19+L32_20.. 
+L32_21+L32_22+L32_23+L32_24+L32_25+L32_26+L32_27+L32_28; 


L33_l=l/(2*L"3)*Blv*m41^2*m42-l/(2*L"3)*Blv*m81^2*m82- 
1/(2*L)*Blv*m51^2*m52*xp2; 

L33_2=-l/(2*L"3)*Blv*m51"2*m52*xp2"3-3/(2*L^2)*Blv*mll"2*ml2*xpl"2; 

L33_3=-l/2*Blv*mll^2*ml2-2/L*Blv*mll"2*ml2*xpl- 

1/(2*L"3)*Blv*mll"2*ml2*xpl^3; 

L33_4=+l/(L^3)*Blv*m41*m81*m82+l/(2*L"3)*Blv*m42*m81"2+l/(L"3)*Blv*m41*m 
81*m82 ; 

L33_5=+l/(2*L"3)*Blv*m42*m81"2-l/(L^3)*Blv*xp2"2*m81*m51*m52; 

L33_6=-l/(2*L^3)*Blv*xp2^2*m82*m51"2+l/(L^3)*Blv*xp2^2*m41*m51*m52; 
L33_7=+l/(2*L"3)*Blv*xp2"2*m42*m51^2-2/(L"2)*Blv*m81*mll*ml2*xpl; 
L33_8=-l/(L^3)*Blv*m81*mll*ml2*xpl"2-l/L*Blv*m81*mll*ml2- 
1/(L^2)*Blv*m82*mll^2*xpl ; 

L33_9=-l/(2*L"3)*Blv*m82*mll"2*xpl"2-l/(2*L)*Blv*m82*mll"2; 

L33_10=-l/(2*L^3)*Blv*xp2*m81^2*m52-l/(L^3)*Blv*xp2*m81*m82*m51; 
L33_ll=-1/(2*L"2)*Blv*m81"2*ml2-l/(2*L"3)*Blv*m81"2*ml2*xpl; 

L33_12=-l/(L"2)*Blv*m81*m82*mll-l/(L"3)*Blv*m81*m82*mll*xpl ; 

L33_13=-l/(2*L^3)*Blv*xp2*m41^2*m52-l/(L^3)*Blv*xp2*m41*m42*m51; 
L33_14=-l/(2*L"2)*Blv*m41"2*ml2-l/(2*L"3)*Blv*m41"2*ml2*xpl- 
1/(L^2)*Blv*m41*m42*mll; 

L33_15=l/(L^3)*Blv*m41*m42*mll*xpl+l/L*Blv*m41*mll*ml2+2/(L"2)*Blv*m41*m 
ll*ml2*xpl; 

L33_16=+l/(L^3)*Blv*m41*mll*ml2*xpl^2+l/(2*L)*Blv*m42*mll^2; 

L33_17=+l/(L"2)*Blv*m42*mll^2*xpl+l/(2*L"3)*Blv*m42*mll^2*xpl"2; 
L33_18=+l/(L^3)*Blv*xp2*m41*m81*m52+l/(L'-B)*Blv*xp2*m41*m82*m51; 
L33_19=+l/(L^3)*Blv*xp2*m42*m81*m51+l/(L"2)*Blv*m41*m81*ml2; 

L33_20=+l/(L^3)*Blv*m41*m81*ml2*xpl+l/(L"2)*Blv*m41*m82*mll; 

L33_21=+l/(L^3)*Blv*m41*m82*mll*xpl+l/(L^2)*Blv*m42*m81*mll; 

L33_22=+l/(L^3)*Blv*m42*m81*mll*xpl-l/(L^3)*Blv*m81*mll*m52*xp2*xpl; 


46 



L33_23=-l/(L^2)*Blv*m81*mll*m52*xp2-l/(L^3)*Blv*m81*ml2*m51*xp2*xpl; 
L33_24 =-1/(L"2)*Blv*m81*ml2*m51*xp2-l/(L^3)*Blv*m82*mll*m51*xp2*xpl; 
L33_25=-l/(L^2)*Blv*m82*mll*m51*xp2+l/(L^3)*Blv*m41*mll*m52*xp2*xpl; 
L33_26=+l/(L^2)*Blv*m41*mll*m52*xp2+l/(L^3)*Blv*m41*ml2*m51*xp2*xpl; 
L33_27=+l/(L"2)*Blv*m41*ml2*m51*xp2+l/(L"3)*Blv*m42*mll*m51*xp2*xpl; 
L33_28=+l/(L^2)*Blv*m42*mll*m51*xp2-l/(L^2)*Blv*mll*m51*m52*xp2"2 ; 
L33_29=-l/(L"3)*Blv*mll*m51*m52*xp2^2*xpl-l/(2*L"2)*Blv*ml2*m51"2*xp2"2 
L33_30=-l/(2*L^3)*Blv*ml2*m51"2*xp2^2*xpl-l/(2*L)*Blv*mll"2*m52*xp2; 
L33_31=-l/(L^2)*Blv*mll^2*m52*xp2*xpl-l/(2*L^3)*Blv*mll^2*m52*xp2*xpl"2 
L33_32=-l/L*Blv*mll*ml2*m51*xp2-2/(L"2)*Blv*mll*ml2*m51*xp2*xpl; 
L33_33=-l/(L^3)*Blv*mll*ml2*m51*xp2*xpl"2- 
Clv/deltasat^2*kpsi^3*mll^2*ml2; 

L33_34=-Clv/deltasat^2*k:v^3*m21^2*m22-Clv/deltasat^2*k:r^3*m31^2*m32 ; 

L33_35=-Clv/deltasat^2*ky"3*m41^2*m42- 

2*Clv/deltasat^2*kpsi^2*kv*m2l*mll*ml2; 

L33_36=-Clv/deltasat^2*kpsi^2*kv*m22*mll^2- 
2*Clv/deltasat^2*kpsi^2*kr*m31*mll*ml2; 

L33_37=-Clv/deltasat^2*kpsi^2*kr*m32*ml1^2- 
2*Clv/deltasat^2*kpsi^2*ky*m4l*mll*ml2; 

L33_38=-Clv/deltasat^2*kpsi^2*ky*m42*mll^2- 
Clv/deltasat^2*kpsi*kv^2*m21^2*ml2 ; 
L33_39=-2*Clv/deltasat"2*kpsi*kv"2*m21*m22*mll- 
Clv/deltasat^2*kpsi*kr^2*m31^2*ml2; 

L33_40=-2*Clv/deltasat^2*kpsi*kr^2*m31*m32*mll- 
Clv/deltasat^2*kpsi*ky^2*m41^2*ml2; 

L33_41=-2*Clv/deltasat^2*kpsi*ky^2*m41*m42*mll- 
2*Clv/deltasat^2*kv^2*kr*m31*m2l*m22; 
L33_42=-Clv/deltasat"2*kv"2*kr*m32*m21"2- 
2*Clv/deltasat^2*kv^2*ky*m4l*m2l*m22; 
L33_43=-Clv/deltasat"2*kv"2*ky*m42*m21"2- 
2*Clv/deltasat^2*kv*kr^2*m2I*m31*m32; 
L33_44=-Clv/deltasat"2*kv*kr"2*m22*m31"2- 
2*Clv/deltasat^2*kv*ky^2*m2l*m4l*m42; 
L33_45=-Clv/deltasat^2*kv*ky^2*m22*m41^2- 
2*Clv/deltasat^2*kr^2*ky*m4I*m31*m32; 

L33_46=-Clv/deltasat^2*kr^2*ky*m42*m31^2- 
2*Clv/deltasat^2*kr*ky^2*m31*m4l*m42; 
L33_47=-Clv/deltasat"2*kr*ky"2*m32*m41"2- 
2*Clv/deltasat^2*kpsi*kv*kr*m31*m2l*ml2; 

L33_48=-2*Clv/deltasat^2*kpsi*kv*kr*m31*m22*ml1; 

L33_4 9=-2*Clv/deltasat^2*kpsi*kv*kr*m32*m2l*ml1; 

L33_50=-2*Clv/deltasat^2*kpsi*kv*ky*m4l*m2l*ml2; 

L33_51=-2*Clv/deltasat^2*kpsi*kv*ky*m4l*m22*ml1; 

L33_52=-2*Clv/deltasat^2*kpsi*kv*ky*m42*m2l*ml1; 

L33_53=-2*Clv/deltasat^2*kpsi*kr*ky*m4I*m31*ml2; 

L33_54=-2*Clv/deltasat^2*kpsi*kr*ky*m4l*m32*ml1; 

L33_55=-2*Clv/deltasat^2*kpsi*kr*ky*m42*m31*ml1; 
L33_56=-2*Clv/deltasat"2*kv*kr*ky*m41*m21*m32- 
2*Clv/deltasat^2*kv*kr*ky*m4I*m22*m31; 
L33_57=-2*Clv/deltasat"2*kv*kr*ky*m42*m21*m31; 


L33 = L33_l+L33_2+L33_3+L33_4+L33_5+L33_6+L33_7+L33_8+L33_9+L33_10... 

+L33_11+L33_12+L33_13+L33_14+L33_15+L33_16+L33_17+L33_18+L33_19+L33_20. 
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+L33_21+L33_22+L33_23+L33_24+L33_25+L33_26+L33_27+L33_28+L33_29+L33_30.. 


+L33_31+L33_32+L33_33+L33_34+L33_35+L33_36+L33_37+L33_38+L33_39+L33_40.. 

+L33_41+L33_42+L33_44+L33_45+L33_46+L33_47+L33_48+L33_49+L33_50+L33_51.. 
+L33_52+L33_53+L33_54+L33_55+L33_56+L33_57; 


L34_l=l/(2*L^3)*Blv*m41*m42^2-l/(2*L^3)*Blv*m81*m82^2- 
1/(2*L)*Blv*m51*m52^2*xp2; 

L34_2=-l/(2*L^3)*Blv*m51*m52^2*xp2^3-3/(2*L^2)*Blv*mll*ml2^2*xpl"2; 

L34_3=-l/2*Blv*mll*ml2^2-2/L*Blv*mll*ml2^2*xpl- 

1/(2*L"3)*Blv*mll*ml2"2*xpl^3; 

L34_4=+l/(2*L^3)*Blv*m41*m82^2+l/(L"3)*Blv*m42*m81*m82+l/(2*L''3)*Blv*m41 
*m82"2; 

L34_5=+l/(L^3)*Blv*m42*m81*m82-l/(2*L^3)*Blv*xp2^2*m81*m52"2; 

L34_6=-l/(L^3)*Blv*xp2"2*m82*m51*m52+l/(2*L^3)*Blv*xp2"2*m41*m52"2; 
L34_7=+l/(L"3)*Blv*xp2"2*m42*m51*m52-l/(L"2)*Blv*m81*ml2"2*xpl; 

L34_8=-l/(2*L^3)*Blv*m81*ml2^2*xpl^2-l/(2*L)*Blv*m81*ml2^2; 

L34_9=-2/(L^2)*Blv*m82*mll*ml2*xpl-l/(L^3)*Blv*m82*mll*ml2*xpl"2; 
L34_10=-l/L*Blv*m82*mll*ml2-l/ (L''3) *Blv*xp2*m81*m82*m52- 
1/(2*L^3)*Blv*xp2*m82^2*m51; 

L34_ll=-1/(L"2)*Blv*m81*m82*ml2-l/(L"3)*Blv*m81*m82*ml2*xpl- 
1/(2*L^2)*Blv*m82^2*mll; 

L34_12=-l/(2*L"3)*Blv*m82"2*mll*xpl-l/(L"3)*Blv*xp2*m41*m42*m52; 
L34_13=-l/(2*L"3)*Blv*xp2*m42^2*m51-l/(L"2)*Blv*m41*m42*ml2; 

L34_14=-l/(L^3)*Blv*m41*m42*ml2*xpl-l/(2*L^2)*Blv*m42^2*mll; 

L34_15=-l/(2*L"3)*Blv*m42"2*mll*xpl+l/(2*L)*Blv*m41*ml2^2+ 

1/(L^2)*Blv*m41*ml2"2*xpl; 

L34_16=+l/(2*L"3)*Blv*m41*ml2"2*xpl^2+l/L*Blv*m42*mll*ml2; 

L34_17=+2/(L"2)*Blv*m42*mll*ml2*xpl+l/(L"3)*Blv*m42*mll*ml2*xpl"2; 
L34_18=+l/(L^3)*Blv*xp2*m41*m82*m52+l/(L^3)*Blv*xp2*m42*m81*m52; 
L34_19=+l/(L"3)*Blv*xp2*m42*m82*m51+l/(L"2)*Blv*m41*m82*ml2; 

L34_20=+l/(L^3)*Blv*m41*m82*ml2*xpl+l/(L^2)*Blv*m42*m81*ml2; 

L34_21=+l/(L"3)*Blv*m42*m81*ml2*xpl+l/(L"2)*Blv*m42*m82*mll; 

L34_22=+l/(L^3)*Blv*m42*m82*mll*xpl-l/(L^3)*Blv*m81*ml2*m52*xp2*xpl; 
L34_23=-l/(L^2)*Blv*m81*ml2*m52*xp2-l/(L^3)*Blv*m82*mll*m52*xp2*xpl; 
L34_24=-l/(L^2)*Blv*m82*mll*m52*xp2-l/(L^3)*Blv*m82*ml2*m51*xp2*xpl; 
L34_25=-l/(L^2)*Blv*m82*ml2*m51*xp2+l/(L^3)*Blv*m41*ml2*m52*xp2*xpl; 
L34_26= +1/(L"2)*Blv*m41*ml2*m52*xp2+l/(L^3)*Blv*m42*mll*m52*xp2*xpl; 
L34_27=+l/(L"2)*Blv*m42*mll*m52*xp2+l/(L"3)*Blv*m42*ml2*m51*xp2*xpl; 
L34_28=+l/(L^2)*Blv*m42*ml2*m51*xp2-l/(2*L"2)*Blv*mll*m52"2*xp2"2; 
L34_29=-l/(2*L^3)*Blv*mll*m52"2*xp2^2*xpl-l/(L"2)*Blv*ml2*m51*m52*xp2"2; 
L34_30=-l/(L^3)*Blv*ml2*m51*m52*xp2"2*xpl-l/L*Blv*mll*ml2*m52*xp2; 
L34_31=-2/(L"2)*Blv*mll*ml2*m52*xp2*xpl- 
1/(L"3)*Blv*mll*ml2*m52*xp2*xpl^2; 

L34_32=-l/(2*L)*Blv*ml2^2*m51*xp2-l/(L"2)*Blv*ml2^2*m51*xp2*xpl; 
L34_33=-l/(2*L^3)*Blv*ml2"2*m51*xp2*xpl^2- 
Clv/deltasat^2*k:psi^3*mll*ml2^2; 

L34_34=-Clv/deltasat^2*k:v^3*m21*m22^2-Clv/deltasat^2*k:r^3*m31*m32^2; 

L34_35=-Clv/deltasat"2*ky"3*m41*m42^2- 

Clv/deltasat^2*kpsi^2*kv*m2l*ml2^2; 


48 



L34_36=-2*Clv/deltasat^2*kpsi^2*kv*m22*mll*ml2- 
Clv/deltasat^2*kpsi^2*kr*m31*ml2^2; 
L34_37=-2*Clv/deltasat"2*kpsi^2*kr*m32*mll*ml2- 
Clv/deltasat^2*kpsi^2*ky*m4l*ml2^2; 
L34_38=-2*Clv/deltasat^2*kpsi^2*ky*m42*mll*ml2; 
L34_39=-2*Clv/deltasat^2*kpsi*kv^2*m21*m22*ml2- 
Clv/deltasat^2*kpsi*kv^2*m22^2*ml1; 
L34_40=-2*Clv/deltasat^2*kpsi*kr^2*m31*m32*ml2- 


Clv/deltasat^2*kpsi* 
L34_4l=-2*Clv/deltas 
Civ/deltasat^2*kpsi* 
L34_42=-Clv/deltasat 
2*Clv/deltasat^2*kv^ 
L34_43=-Clv/deltasat 
2*Clv/deltasat^2*kv 
L34_44=-Clv/deltasat 
2*Clv/deltasat^2*kv 
L34_45 = -Clv/deltas 
2*Clv/deltasat^2*kv* 
L34_46 = -Clv/deltas 
2*Clv/deltasat^2*kr 
L34_47 = -Clv/deltas 
2*Clv/deltasat^2*kr* 


L34_48 = -2*Clv/delt 
L34_49 = -2*Clv/delt 
L34_50 = -2*Clv/delt 
L34_51 = -2*Clv/delt 
L34_52 = -2*Clv/delt 
L34_53 = -2*Clv/delt 
L34_54 = -2*Clv/delt 
L34_55 = -2*Clv/delt 
L34_56 = -2*Clv/delt 
L34_57 = -2*Clv/delt 
2*Clv/deltasat^2*kv* 
L34_58 = -2*Clv/delt 


kr^2*m32^2*mll; 

at^2*kpsi*ky^2*m41*m42*ml2- 

ky^2*m42^2*mll; 

"2*kv^2*kr*m31*m22"2- 

2*kr*m32*m21*m22; 

2*kv"2*ky*m41*m22^2- 

2*ky*m42*m21*m22; 

2*kv*kr^2*m21*m32^2- 
kr"2*m22*m31*m32; 
at^2*kv*ky^2*m21*m42^2- 
ky^2*m22*m41*m42; 
at"2*kr^2*ky*m41*m32"2- 
2*ky*m42*m31*m32; 
at"2*kr*ky"2*m31*m42"2- 
ky^2*m32*m41*m42; 
asat^2*kpsi*kv*kr*m31*m22*ml2; 
asat^2*kpsi*kv*kr*m32*m21*ml2; 
asat^2*kpsi*kv*kr*m32*m22*mll; 
asat^2*kpsi*kv*ky*m41*m22*ml2; 
asat^2*kpsi*kv*ky*m42*m21*ml2; 
asat^2*kpsi*kv*ky*m42*m22*mll; 
asat^2*kpsi*kr*ky*m41*m32*ml2; 
asat^2*kpsi*kr*ky*m42*m31*ml2; 
asat^2*kpsi*kr*ky*m42*m32*mll; 
asat^2*kv*kr*ky*m41*m22*m32- 
kr*ky*m42*m21*m32; 
asat^2*kv*kr*ky*m42*m22*m31 ; 


L34 = L34_l+L34_2+L34_3+L34_4+L34_5+L34_6+L34_7+L34_8+L34_9+L34_10.. 
+L34_11+L34_12+L34_13+L34_14+L34_15+L34_16+L34_17+L34_18+L34_19+L34_20 
+L34_21+L34_22+L34_23+L34_24+L34_25+L34_26+L34_27+L34_28+L34_29+L34_30 
+L34_31+L34_32+L34_33+L34_34+L34_35+L34_36+L34_37+L34_38+L34_39+L34_40 
+L34_41+L34_42+L34_44+L34_45+L34_46+L34_47+L34_48+L34_49+L34_50+L34_51 
+L34_52+L34_53+L34_54+L34_55+L34_56+L34_57+L34_58; 


L41_l=l/(6*L"'3)*Blr*m41"3-l/(6*L^3)*Blr*m81"3-l/(6*L)*Blr*m51"3*xp2; 
L41_2=-l/(6*L^3)*Blr*m51^3*xp2"3-l/(2*L"2)*Blr*mll^3*xpl^2- 
l/6*Blr*mll^3; 

L41_3=-2/(3*L)*Blr*mll^3*xpl- 

1/(6*L"3)*Blr*mll"3*xpl^3+l/(2*L"3)*Blr*m41*m81^2; 
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L41_4=+l/(2*L^3)*Blr*m41*m81^2-l/(2*L^3)*Blr*xp2^2*m81*m51^2; 

L41_5 = +1/(2*L^3)*Blr*xp2"2*m41*m51"2-l/(L^2)*Blr*m81*mll"2*xpl; 

L41_6 = -1/(2*L^3)*Blr*m81*mll"2*xpl^2-l/(2*L)*Blr*m81*mll^2; 

L41_7 = -1/(2*L^3)*Blr*xp2*m81"2*m51-l/(2*L^2)*Blr*m81"2*mll; 

L41_8 = -1/(2*L^3)*Blr*m81"2*mll*xpl-l/(2*L^3)*Blr*xp2*m41"2*m51; 

L41_9 = -1/(2*L"2)*Blr*m41^2*mll- 

1/(2*L"3)*Blr*m41"2*mll*xpl+l/(2*L)*Blr*m41*mll^2; 

L41_10 = +1/(L"2)*Blr*m41*mll"2*xpl+l/(2*L^3)*Blr*m41*mll^2*xpl"2; 

L41_ll = +1/(L'-B)*Blr*xp2*m41*m81*m51+l/(L^2)*Blr*m41*m81*mll; 

L41_12 = +1/(L^3)*Blr*m41*m81*mll*xpl-l/(L"3)*Blr*m81*mll*m51*xp2*xpl; 
L41_13 = -1/(L^2)*Blr*m81*mll*m51*xp2+l/(L^3)*Blr*m41*mll*m51*xp2*xpl; 
L41_14 = +1/(L^2)*Blr*m41*mll*m51*xp2-l/(2*L"2)*Blr*mll*m51"2*xp2"2; 
L41_15 = -1/(2*L^3)*Blr*mll*m51^2*xp2^2*xpl-l/(2*L)*Blr*mll"2*m51*xp2; 
L41_16 = -1/(L"2)*Blr*mll"2*m51*xp2*xpl- 
1/(2*L^3)*Blr*mll"2*m51*xp2*xpl"2; 

L41_17 = -l/3*Clr/deltasat^2*kpsi^3*mll"3-l/3*Clr/deltasat^2*kv"3*m21"3; 
L41_18 = -l/3*Clr/deltasat"2*kr^3*m31"3-l/3*Clr/deltasat"2*ky"3*m41"3; 
L41_19 = -Clr/deltasat^2*kpsi^2*kv*m21*mll^2- 
Clr/deltasat^2*kpsi^2*kr*m31*ml1^2; 

L41_20 = -Clr/deltasat^2*kpsi^2*ky*m41*mll^2- 
Clr/deltasat^2*kpsi*kv^2*m21^2*ml1; 

L41_21 = -Clr/deltasat^2*kpsi*kr^2*m31^2*mll- 
Clr/deltasat^2*kpsi*ky^2*m41^2*ml1; 

L41_22 = -Clr/deltasat"2*kv^2*kr*m31*m21^2- 
Clr/deltasat^2*kv^2*ky*m4l*m21^2; 

L41_23 = -Clr/deltasat^2*kv*kr"'2*m21*m31^2- 
Clr/deltasat^2*kv*ky^2*m2l*m41^2; 

L41_24 = -Clr/deltasat"2*kr^2*ky*m41*m31"2- 
Clr/deltasat^2*kr*ky^2*m31*m41^2; 

L41_25 = -2*Clr/deltasat^2*kpsi*kv*kr*m31*m2l*ml1; 

L41_26 = -2*Clr/deltasat^2*kpsi*kv*ky*m4l*m2l*ml1; 

L41_27 = -2*Clr/deltasat^2*kpsi*kr*ky*m4l*m31*ml1; 

L41_28 = -2*Clr/deltasat"2*kv*kr*ky*m41*m21*m31; 


L41=L41_1+L41_2+L41_3+L41_4+L41_5+L41_6+L41_7+L41_8+L41_9+L41_10... 

+L41_11+L41_12+L41_13+L41_14+L41_15+L41_16+L41_17+L41_18+L41_19+L41_20.. 

+L41_21+L41_22+L41_23+L41_24+L41_25+L41_26+L41_27+L41_28; 

L42_l=l/(6*L"3)*Blr*m42^3-l/(6*L^3)*Blr*m82^3-l/(6*L)*Blr*m52"3*xp2 ; 
L42_2=-l/(6*L^3)*Blr*m52^3*xp2^3-l/(2*L^2)*Blr*ml2^3*xpl^2-l/6*Blr*ml2"3 
L42_3=-2/(3*L)*Blr*ml2^3*xpl- 

1/(6*L^3)*Blr*ml2^3*xpl"3+l/(2*L"3)*Blr*m42*m82"2; 

L42_4=+l/(2*L"3)*Blr*m42*m82"2-l/(2*L"3)*Blr*xp2"2*m82*m52"2 ; 

L42_5=+l/(2*L"3)*Blr*xp2"2*m42*m52"2-l/(L"2)*Blr*m82*ml2"2*xpl; 

L42_6=-l/(2*L^3)*Blr*m82*ml2^2*xpl^2-l/(2*L)*Blr*m82*ml2^2 ; 

L42_7=-l/(2*L"3)*Blr*xp2*m82"2*m52-l/(2*L"2)*Blr*m82"2*ml2; 

L42_8=-l/(2*L"3)*Blr*m82^2*ml2*xpl-l/(2*L^3)*Blr*xp2*m42^2*m52; 

L42_9=-l/(2*L^2)*Blr*m42^2*ml2- 

1/(2*L^3)*Blr*m42"2*ml2*xpl+l/(2*L)*Blr*m42*ml2^2; 

L42_10=+l/(L^2)*Blr*m42*ml2^2*xpl+l/(2*L^3)*Blr*m42*ml2^2*xpl"2 ; 

L42_ll = +1/(L"3)*Blr*xp2*m42*m82*m52+l/(L"2)*Blr*m42*m82*ml2; 
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L42_12 = +1/(L"3)*Blr*m42*m82*ml2*xpl-l/(L"3)*Blr*m82*ml2*m52*xp2*xpl; 

L42_13 = -1/(L"2)*Blr*m82*ml2*m52*xp2+l/(L"3)*Blr*m42*ml2*m52*xp2*xpl; 
L42_14=+l/(L^2)*Blr*m42*ml2*m52*xp2-l/(2*L^2)*Blr*ml2*m52"2*xp2"2 ; 

L42_15=-l/(2*L^3)*Blr*ml2*m52^2*xp2^2*xpl-l/(2*L)*Blr*ml2"2*m52*xp2; 
L42_16=-l/(L^2)*Blr*ml2^2*m52*xp2*xpl-l/(2*L"3)*Blr*ml2^2*m52*xp2*xpl"2; 
L42_17=-l/3*Clr/deltasat^2*kpsi"3*ml2^3-l/3*Clr/deltasat^2*kv"3*m22"3; 
L42_18 =-l/3*Clr/deltasat"2*kr"3*m32"3-l/3*Clr/deltasat^2*ky"3*m42"3 ; 

L42_19 = -Clr/deltasat^2*kpsi^2*kv*m22*ml2^2- 
Clr/deltasat^2*kpsi^2*kr*m32*ml2^2; 

L42_20 = -Clr/deltasat^2*kpsi^2*ky*m42*ml2^2- 
Clr/deltasat^2*kpsi*kv^2*m22^2*ml2; 

L42_21 = -Clr/deltasat"2*kpsi*kr"2*m32"2*ml2- 
Clr/deltasat^2*kpsi*ky^2*m42^2*ml2; 

L42_22 = -Clr/deltasat^2*kv^2*kr*m32*m22^2- 
Clr/deltasat^2*kv^2*ky*m42*m22^2 ; 

L42_23 = -Clr/deltasat^2*kv*kr^2*m22*m32^2- 
Clr/deltasat^2*kv*ky^2*m22*m42^2; 

L42_24 =-Clr/deltasat"2*kr"2*ky*m42*m32^2- 
Clr/deltasat^2*kr*ky^2*m32*m42^2; 

L42_25 = -2*Clr/deltasat^2*kpsi*kv*kr*m32*m22*ml2; 

L42_26 = -2*Clr/deltasat^2*kpsi*kv*ky*m42*m22*ml2; 

L42_27 = -2*Clr/deltasat^2*kpsi*kr*ky*m42*m32*ml2; 

L42_28 = -2*Clr/deltasat^2*kv*kr*ky*m42*m22*m32; 


L42=L42_1+L42_2+L42_3+L42_4+L42_5+L42_6+L42_7+L42_8+L42_9+L42_10... 
+L42_11+L42_12+L42_13+L42_14+L42_15+L42_16+L42_17+L42_18+L42_19+L42_20.. 
+L42_21+L42_22+L42_23+L42_24+L42_25+L42_26+L42_27+L42_28; 


L43_l=l/(2*L"3)*Blr*m41^2*m42-l/(2*L"3)*Blr*m81^2*m82- 
1/(2*L)*Blr*m51^2*m52*xp2; 

L43_2=-l/(2*L"3)*Blr*m51"2*m52*xp2^3-3/(2*L^2)*Blr*mll"2*ml2*xpl"2 ; 

L43_3=-l/2*Blr*mll^2*ml2-2/L*Blr*mll"2*ml2*xpl- 

1/(2*L^3)*Blr*mll"2*ml2*xpl^3; 

L43_4=+l/(L^3)*Blr*m41*m81*m82+l/(2*L"3)*Blr*m42*m81"2+l/(L"3)*Blr*m41*m 
81*m82 ; 

L43_5=+l/(2*L"3)*Blr*m42*m81"2-l/(L^3)*Blr*xp2"2*m81*m51*m52; 

L43_6=-l/(2*L^3)*Blr*xp2^2*m82*m51^2+l/(L^3)*Blr*xp2^2*m41*m51*m52 ; 

L43_7=+l/(2*L"3)*Blr*xp2"2*m42*m51"2-2/(L"2)*Blr*m81*mll*ml2*xpl; 
L43_8=-l/(L^3)*Blr*m81*mll*ml2*xpl"2-l/L*Blr*m81*mll*ml2- 
1/(L^2)*Blr*m82*mll^2*xpl ; 

L43_9=-l/(2*L"3)*Blr*m82*mll"2*xpl^2-l/(2*L)*Blr*m82*mll"2; 

L43_10=-l/(2*L^3)*Blr*xp2*m81^2*m52-l/(L"3)*Blr*xp2*m81*m82*m51 ; 

L43_ll=-1/(2*L^2)*Blr*m81"2*ml2-l/(2*L"3)*Blr*m81"2*ml2*xpl- 
1/(L"2)*Blr*m81*m82*mll; 

L43_12=-l/(L^3)*Blr*m81*m82*mll*xpl-l/(2*L^3)*Blr*xp2*m41^2*m52 ; 

L43_13=-l/(L^3)*Blr*xp2*m41*m42*m51-l/(2*L"2)*Blr*m41"2*ml2; 

L43_14=-l/(2*L^3)*Blr*m41^2*ml2*xpl-l/(L"2)*Blr*m41*m42*mll ; 

L43_15=l/(L"3)*Blr*m41*m42*mll*xpl+l/L*Blr*m41*mll*ml2+2/(L"2)*Blr*m41*m 
ll*ml2*xpl; 

L43_16=+l/(L^3)*Blr*m41*mll*ml2*xpl^2+l/(2*L)*Blr*m42*mll"2 ; 
L43_17=+l/(L^2)*Blr*m42*mll^2*xpl+l/(2*L"3)*Blr*m42*mll^2*xpl"2; 
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L43_18=+l/(L"3)*Blr*xp2*m41*m81*m52+l/(L^3)*Blr*xp2*m41*m82*m51 ; 
L43_19=+l/(L"3)*Blr*xp2*m42*m81*m51+l/(L"2)*Blr*m41*m81*ml2; 

L43_20=+l/(L^3)*Blr*m41*m81*ml2*xpl+l/(L^2)*Blr*m41*m82*mll ; 

L43_21=+l/(L^3)*Blr*m41*m82*mll*xpl+l/(L^2)*Blr*m42*m81*mll; 

L43_22=+l/(L"3)*Blr*m42*m81*mll*xpl-l/(L"3)*Blr*m81*mll*m52*xp2*xpl ; 

L43_23=-l/(L^2)*Blr*m81*mll*m52*xp2-l/(L^3)*Blr*m81*ml2*m51*xp2*xpl; 
L43_24=-l/(L"2)*Blr*m81*ml2*m51*xp2-l/(L"3)*Blr*m82*mll*m51*xp2*xpl ; 
L43_25=-l/(L"2)*Blr*m82*mll*m51*xp2+l/(L^3)*Blr*m41*mll*m52*xp2*xpl; 
L43_26=+l/(L"2)*Blr*m41*mll*m52*xp2+l/(L^3)*Blr*m41*ml2*m51*xp2*xpl; 
L43_27=+l/(L^2)*Blr*m41*ml2*m51*xp2+l/(L"3)*Blr*m42*mll*m51*xp2*xpl; 
L43_28=+l/(L^2)*Blr*m42*mll*m51*xp2-l/(L"2)*Blr*mll*m51*m52*xp2"2 ; 

L43_29=-l/(L^3)*Blr*mll*m51*m52*xp2^2*xpl-l/(2*L^2)*Blr*ml2*m51"2*xp2"2; 
L43_30=-l/(2*L^3)*Blr*ml2*m51^2*xp2^2*xpl-l/(2*L)*Blr*mll^2*m52*xp2 ; 

L43_31=-l/(L^2)*Blr*mll^2*m52*xp2*xpl-l/(2*L^3)*Blr*mll^2*m52*xp2*xpl"2; 
L43_32=-l/L*Blr*mll*ml2*m51*xp2-2/(L"2)*Blr*mll*ml2*m51*xp2*xpl ; 

L43_33=-l/(L^3)*Blr*mll*ml2*m51*xp2*xpl"2- 
Clr/deltasat^2*kpsi^3*mll^2*ml2; 

L4 3_34=-Clr/deltasat^2*kv^3*m21^2*m22-Clr/deltasat^2*kr^3*m31^2*m32; 

L4 3_35=-Clr/deltasat"'2*ky^3*m41^2*m42- 
2*Clr/deltasat^2*kpsi^2 *kv*m2l*mll*ml2; 

L43_36=-Clr/deltasat^2*kpsi^2*kv*m22*mll^2- 
2*Clr/deltasat^2*kpsi^2*kr*m31*mll*ml2 ; 

L43_37=-Clr/deltasat^2*kpsi^2*kr*m32*mll^2- 
2*Clr/deltasat^2*kpsi^2*ky*m4l*mll*ml2; 

L43_38=-Clr/deltasat^2*kpsi^2*ky*m42*mll^2- 
Clr/deltasat^2*kpsi*kv^2*m21^2*ml2 ; 

L43_39 = -2*Clr/deltasat^2*kpsi*kv^2*m21*m22*mll- 
Clr/deltasat^2*kpsi*kr^2*m31^2*ml2; 

L43_40 = -2*Clr/deltasat^2*kpsi*kr^2*m31*m32*mll- 
Clr/deltasat^2*kpsi*ky^2*m41^2*ml2 ; 

L43_41 =-2*Clr/deltasat^2*kpsi*ky^2*m41*m42*mll- 
2*Clr/deltasat^2*kv^2*kr*m31*m21*m22 ; 

L43_42 = -Clr/deltasat"2*kv^2*kr*m32*m21"2- 
2*Clr/deltasat^2*kv^2*ky*m41*m21*m22 ; 

L43_43 = -Clr/deltasat"2*kv^2*ky*m42*m21"2- 
2*Clr/deltasat^2*kv*kr^2*m2I*m31*m32; 

L43_44 = -Clr/deltasat^2*kv*kr^2*m22*m31"2- 
2*Clr/deltasat^2*kv*ky^2*m2l*m4l*m42; 

L43_45 = -Clr/deltasat^2*kv*ky^2*m22*m41^2- 
2*Clr/deltasat^2*kr^2*ky*m4I*m31*m32; 

L43_46 = -Clr/deltasat"2*kr^2*ky*m42*m31^2- 
2*Clr/deltasat^2*kr*ky^2*m31*m4l*m42; 

L43_47 = -Clr/deltasat^2*kr*ky^2*m32*m41^2- 
2*Clr/deltasat^2*kpsi*kv*kr*m31*m2l*ml2; 

L43_48 = -2*Clr/deltasat^2*kpsi*kv*kr*m31*m22*mll ; 

L43_49 = -2*Clr/deltasat^2*kpsi*kv*kr*m32*m2l*ml1; 

L43_50 = -2*Clr/deltasat^2*kpsi*kv*ky*m41*m21*ml2 ; 

L43_51 = -2*Clr/deltasat^2*kpsi*kv*ky*m4l*m22*ml1; 

L43_52 = -2*Clr/deltasat^2*kpsi*kv*ky*m42*m2l*ml1; 

L43_53 = -2*Clr/deltasat^2*kpsi*kr*ky*m4I*m31*ml2; 

L43_54 = -2*Clr/deltasat^2*kpsi*kr*ky*m41*m32*mll ; 

L43_55 = -2*Clr/deltasat^2*kpsi*kr*ky*m42*m31*ml1; 

L43_56 = -2*Clr/deltasat"2*kv*kr*ky*m41*m21*m32- 
2*Clr/deltasat^2*kv*kr*ky*m41*m22*m31 ; 

L43_57 = -2*Clr/deltasat"2*kv*kr*ky*m42*m21*m31; 
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L43= L43_l+L43_2+L43_3+L43_4+L43_5+L43_6+L43_7+L43_8+L43_9+L43_10... 

+L43_11+L43_12+L43_13+L43_14+L43_15+L43_16+L43_17+L43_18+L43_19+L43_20.. 

+L43_21+L43_22+L43_23+L43_24+L43_25+L43_26+L43_27+L43_28+L43_29+L43_30.. 

+L43_31+L43_32+L43_33+L43_34+L43_35+L43_36+L43_37+L43_38+L43_39+L43_40.. 

+L43_41+L43_42+L43_44+L43_45+L43_46+L43_47+L43_48+L43_49+L43_50+L43_51.. 
L43_52+L43_53+L43_54+L43_55+L43_56+L43_57; 


L44_l=l/(2*L"3)*Blr*m41*m42^2-l/(2*L"3)*Blr*m81*m82^2- 
1/(2*L)*Blr*m51*m52"2*xp2; 

L44_2=-l/(2*L"3)*Blr*m51*m52"2*xp2"3-3/(2*L^2)*Blr*mll*ml2"2*xpl"2; 

L44_3=-l/2*Blr*mll*ml2"2-2/L*Blr*mll*ml2^2*xpl- 

1/(2*L^3)*Blr*mll*ml2^2*xpl^3; 

L44_4=+l/(2*L"3)*Blr*m41*m82"2+1/(L"3)*Blr*m42*m81*m82+l/(2*L"3)*Blr*m41 
*m82"2; 

L44_5=+l/(L"3)*Blr*m42*m81*m82-l/(2*L"3)*Blr*xp2"2*m81*m52"2; 

L44_6=-l/(L"3)*Blr*xp2"2*m82*m51*m52+l/(2*L"3)*Blr*xp2"2*m41*m52"2; 
L44_7=+l/(L"3)*Blr*xp2"2*m42*m51*m52-l/(L"2)*Blr*m81*ml2"2*xpl; 

L44_8=-l/(2*L"3)*Blr*m81*ml2"2*xpl"2-l/(2*L)*Blr*m81*ml2"2; 

L44_9=-2/(L"2)*Blr*m82*mll*ml2*xpl-l/(L"3)*Blr*m82*mll*ml2*xpl"2; 
L44_10=-l/L*Blr*m82*mll*ml2-l/(L"3)*Blr*xp2*m81*m82*m52- 
1/(2*L"3)*Blr*xp2*m82"2*m51; 

L44_ll=-1/(L"2)*Blr*m81*m82*ml2-l/(L"3)*Blr*m81*m82*ml2*xpl- 
1/(2*L"2)*Blr*m82"2*mll ; 

L44_12=-l/(2*L"3)*Blr*m82"2*mll*xpl-l/(L"3)*Blr*xp2*m41*m42*m52; 
L44_13=-l/(2*L"3)*Blr*xp2*m42"2*m51-l/(L"2)*Blr*m41*m42*ml2; 

L44_14=-l/(L"3)*Blr*m41*m42*ml2*xpl-l/(2*L"2)*Blr*m42"2*mll; 

L44_15=-l/(2*L"3)*Blr*m42"2*mll*xpl+l/(2*L)*Blr*m41*ml2"2+ 

1/(L"2)*Blr*m41*ml2"2*xpl; 

L44_16=+l/(2*L"3)*Blr*m41*ml2"2*xpl"2+l/L*Blr*m42*mll*ml2; 

L44_17=+2/(L"2)*Blr*m42*ml1*ml2*xpl+l/(L"3)*Blr*m42*mll*ml2*xpl"2; 
L44_18=+l/(L"3)*Blr*xp2*m41*m82*m52+l/(L"3)*Blr*xp2*m42*m81*m52; 
L44_19=+l/(L"3)*Blr*xp2*m42*m82*m51+l/(L"2)*Blr*m41*m82*ml2; 

L44_20=+l/(L"3)*Blr*m41*m82*ml2*xpl+l/(L"2)*Blr*m42*m81*ml2; 

L44_21=+l/(L"3)*Blr*m42*m81*ml2*xpl+l/(L"2)*Blr*m42*m82*mll; 

L44_22=+l/(L"3)*Blr*m42*m82*mll*xpl-l/(L"3)*Blr*m81*ml2*m52*xp2*xpl; 
L44_23=-l/(L"2)*Blr*m81*ml2*m52*xp2-l/(L"3)*Blr*m82*mll*m52*xp2*xpl; 
L44_24=-l/(L"2)*Blr*m82*ml1*m52*xp2-l/(L"3)*Blr*m82*ml2*m51*xp2*xpl; 
L44_25=-l/(L"2)*Blr*m82*ml2*m51*xp2+l/(L"3)*Blr*m41*ml2*m52*xp2*xpl; 
L44_26=+l/(L"2)*Blr*m41*ml2*m52*xp2+l/(L"3)*Blr*m42*mll*m52*xp2*xpl; 
L44_27=+l/(L"2)*Blr*m42*ml1*m52*xp2+l/(L"3)*Blr*m42*ml2*m51*xp2*xpl; 
L44_28=+l/(L"2)*Blr*m42*ml2*m51*xp2-l/(2*L"2)*Blr*mll*m52"2*xp2"2; 
L44_29=-l/(2*L"3)*Blr*mll*m52"2*xp2"2*xpl-l/(L"2)*Blr*ml2*m51*m52*xp2"2; 
L44_30=-l/(L"3)*Blr*ml2*m51*m52*xp2"2*xpl-l/L*Blr*mll*ml2*m52*xp2; 
L44_31=-2/(L"2)*Blr*mll*ml2*m52*xp2*xpl- 
1/(L"3)*Blr*mll*ml2*m52*xp2*xpl"2; 

L44_32=-l/(2*L)*Blr*ml2"2*m51*xp2-l/(L"2)*Blr*ml2"2*m51*xp2*xpl; 
L44_33=-l/(2*L"3)*Blr*ml2"2*m51*xp2*xpl"2- 
Clr/deltasat"2*kpsi"3*mll*ml2"2; 
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L44_34 =-Clr/deltasat^2*kv^3*m21*m22^2-Clr/deltasat^2*kr^3*m31*m32^2; 

L44_35=-Clr/deltasat"2*ky"3*m41*m42^2- 

Clr/deltasat^2*kpsi^2*kv*m2l*ml2^2; 

L44_36=-2*Clr/deltasat^2*kpsi^2*kv*m22*mll*ml2- 
Clr/deltasat^2*kpsi^2*kr*m31*ml2^2; 
L44_37=-2*Clr/deltasat^2*kpsi^2*kr*m32*mll*ml2- 
Clr/deltasat^2*kpsi^2*ky*m4l*ml2^2; 

L44_38=-2*Clr/deltasat^2*kpsi^2*ky*m42*mll*ml2; 

L44_39=-2*Clr/deltasat^2*kpsi*kv^2*m21*m22*ml2- 
Clr/deltasat^2*kpsi*kv^2*m22^2*ml1; 

L44_40 = -2*Clr/deltasat^2*kpsi*kr^2*m31*m32*ml2- 
Clr/deltasat^2*kpsi*kr^2*m32^2*ml1; 

L44_41=-2*Clr/deltasat^2*kpsi*ky^2*m41*m42*ml2- 
Clr/deltasat^2*kpsi*ky^2*m42^2*ml1; 

L44_42=-Clr/deltasat^2*kv^2*kr*m31*m22^2- 
2*Clr/deltasat^2*kv^2*kr*m32*m2l*m22; 

L44_43=-Clr/deltasat^2*kv^2*ky*m41*m22^2- 
2*Clr/deltasat^2*kv^2*ky*m42*m2l*m22; 
L44_44=-Clr/deltasat^2*kv*kr^2*m21*m32"2- 
2*Clr/deltasat^2*kv*kr^2*m22*m31*m32; 
L44_45=-Clr/deltasat^2*kv*ky^2*m21*m42^2- 
2*Clr/deltasat^2*kv*ky^2*m22*m4l*m42; 

L44_46=-Clr/deltasat^2*kr^2*ky*m41*m32^2- 
2*Clr/deltasat^2*kr^2*ky*m42*m31*m32; 
L44_47=-Clr/deltasat"2*kr*ky"2*m31*m42^2- 
2*Clr/deltasat^2*kr*ky^2*m32*m4l*m42; 

L44_48=-2*Clr/deltasat^2*kpsi*kv*kr*m31*m22*ml2; 

L44_49 = -2*Clr/deltasat^2*kpsi*kv*kr*m32*m2l*ml2; 

L44_50 =-2*Clr/deltasat^2*kpsi*kv*kr*m32*m22*ml1; 

L44_51 = -2*Clr/deltasat^2*kpsi*kv*ky*m4I*m22*ml2; 

L44_52 = -2*Clr/deltasat^2*kpsi*kv*ky*m42*m2l*ml2; 

L44_53 = -2*Clr/deltasat^2*kpsi*kv*ky*m42*m22*ml1; 

L44_54 =-2*Clr/deltasat^2*kpsi*kr*ky*m4I*m32*ml2; 

L44_55 = -2*Clr/deltasat^2*kpsi*kr*ky*m42*m31*ml2; 

L44_56 = -2*Clr/deltasat^2*kpsi*kr*ky*m42*m32*ml1; 

L44_57 = -2*Clr/deltasat^2*kv*kr*ky*m41*m22*m32- 
2*Clr/deltasat^2*kv*kr*ky*m42*m2l*m32; 

L44_58 =-2*Clr/deltasat^2*kv*kr*ky*m42*m22*m31; 


L44 = L44_l+L44_2+L44_3+L44_4+L44_5+L44_6+L44_7+L44_8+L44_9+L44_10.. 
+L44_11+L44_12+L44_13+L44_14+L44_15+L44_16+L44_17+L44_18+L44_19+L44_20 
+L44_21+L44_22+L44_23+L44_24+L44_25+L44_26+L44_27+L44_28+L44_29+L44_30 
+L44_31+L44_32+L44_33+L44_34+L44_35+L44_36+L44_37+L44_38+L44_39+L44_40 
+L44_41+L44_42+L44_44+L44_45+L44_46+L44_47+L44_48+L44_49+L44_50+L44_51 
+L44_52+L44_53+L44_54+L44_55+L44_56+L44_57+L44_58; 
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L51=l/6*mll^3-l/6*m51"3-l/2*m61*m51^21/2*m21*mll"2; 

L52=l/6*ml2^3-l/6*m52^3-l/2*m62*m52^2+l/2*m22*ml2"2; 

L53=l/2*mll^2*ml2-l/2*m51^2*m52-m61*m51*m52+m21*mll*ml2- 

l/2*m62*m51^2+l/2*m22*mll"2; 

L54=l/2*mll*ml2"2-l/2*m51*m52^2-l/2*m61*m52"2- 

m62*m51*m52+l/2*m21*ml2"2+m22*mll*ml2; 


Rll 
R13 
R22 
R2 4 


N11*L11 + N12*L21 
N11*L14 + N12*L24 
N21*L13 + N22*L23 
N21*L12 + N22*L22 


+ N13*L31 + N14*L41 + N15*L51; 
+ N13*L34 + N14*L44 + N15*L54; 
+ N23*L33 + N24*L43 + N25*L53; 
+ N23*L32 + N24*L42 + N25*L52; 


Kcoef(index)=(1/8)*(3*R11+R13+R22+3*R24); 


end 


end 


end 


end 


% Plot (L, crit TC) curve 
"6 

plot(L_v,Kcoefb'),xlabel('L'),ylabel('{K}'),title('Constant T 
%end 


),grid 
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